Есть ли интеллектуальная и пространственно-эффективная симметричная матрица в numpy, которая автоматически (и прозрачно) заполняет позицию в [j][i], когда [i][j] записывается на?
import numpy
a = numpy.symmetric((3, 3))
a[0][1] = 1
a[1][0] == a[0][1]
# True
print(a)
# [[0 1 0], [1 0 0], [0 0 0]]
assert numpy.all(a == a.T) # for any symmetric matrix
Автоматический эрмитов тоже будет приятным, хотя мне и не понадобится это на момент написания.
Если вы можете позволить себе симметризировать матрицу непосредственно перед выполнением вычислений, следующее должно быть достаточно быстрым:
def symmetrize(a):
return a + a.T - numpy.diag(a.diagonal())
Это работает при разумных предположениях (например, не делая как a[0, 1] = 42, так и противоречивый a[1, 0] = 123 перед запуском symmetrize).
Если вам действительно нужна прозрачная симметризация, вы можете рассмотреть подклассификацию numpy.ndarray и просто переопределить __setitem__:
class SymNDArray(numpy.ndarray):
def __setitem__(self, (i, j), value):
super(SymNDArray, self).__setitem__((i, j), value)
super(SymNDArray, self).__setitem__((j, i), value)
def symarray(input_array):
"""
Returns a symmetrized version of the array-like input_array.
Further assignments to the array are automatically symmetrized.
"""
return symmetrize(numpy.asarray(input_array)).view(SymNDArray)
# Example:
a = symarray(numpy.zeros((3, 3)))
a[0, 1] = 42
print a # a[1, 0] == 42 too!
(или эквивалент с матрицами вместо массивов, в зависимости от ваших потребностей). Этот подход даже обрабатывает более сложные назначения, такие как a[:, 1] = -1, который правильно устанавливает элементы a[1, :].
Обратите внимание, что Python 3 удалил возможность записи def …(…, (i, j),…), поэтому перед запуском с Python 3 код
Более общий вопрос об оптимальном лечении симметричных матриц в numpy тоже прослушивал меня.
Посмотрев на него, я думаю, что ответ, вероятно, заключается в том, что numpy несколько ограничивается поддержкой макета памяти под базовыми процедурами BLAS для симметричных матриц.
Хотя некоторые процедуры BLAS используют симметрию для ускорения вычислений на симметричных матрицах, они по-прежнему используют ту же структуру памяти, что и полная матрица, то есть n^2 space, а не n(n+1)/2. Просто им говорят, что матрица симметрична и использует только значения в верхнем или нижнем треугольнике.
Некоторые из подпрограмм scipy.linalg принимают флаги (например, sym_pos=True on linalg.solve), которые передаются в процедуры BLAS, хотя большая поддержка этого в numpy будет приятной, в частности, обертки для подпрограмм, таких как DSYRK ( симметричное ранговое обновление k), что позволило бы вычислить матрицу Грама быстрее, чем точка (MT, M).
(Может показаться странным беспокоиться об оптимизации для 2x постоянного фактора во времени и/или пространстве, но это может повлиять на этот порог того, насколько большой проблемой вы можете справиться на одной машине...)
Существует ряд хорошо известных способов хранения симметричных матриц, поэтому им не нужно занимать n ^ 2 элемента хранения. Кроме того, можно переписать общие операции для доступа к этим пересмотренным средствам хранения. Окончательной работой являются Голуб и Ван Лоан, Matrix Computations, 3-е издание 1996 года, Университетский университет Джона Хопкинса, разделы 1.27-1.2.9. Например, цитируя их из формы (1.2.2), в симметричной матрице нужно хранить A = [a_{i,j} ] для i >= j. Тогда, предполагая, что вектор, удерживающий матрицу, обозначается V, а A - n-by-n, поместите a_{i,j} в
V[(j-1)n - j(j-1)/2 + i]
Это предполагает 1-индексирование.
Голуб и Ван Лоан предлагают Алгоритм 1.2.3, который показывает, как получить доступ к такой сохраненной V, чтобы вычислить y = V x + y.
Голуб и Ван Лоан также обеспечивают способ хранения матрицы в диагональной доминирующей форме. Это не сохраняет память, но поддерживает доступ к ним для некоторых других видов операций.
Это простой python, но не numpy, но я просто скомпилировал рутину для заполнения симметричную матрицу (и тестовую программу, чтобы убедиться, что она правильная):
import random
# fill a symmetric matrix with costs (i.e. m[x][y] == m[y][x]
# For demonstration purposes, this routine connect each node to all the others
# Since a matrix stores the costs, numbers are used to represent the nodes
# so the row and column indices can represent nodes
def fillCostMatrix(dim): # square array of arrays
# Create zero matrix
new_square = [[0 for row in range(dim)] for col in range(dim)]
# fill in main diagonal
for v in range(0,dim):
new_square[v][v] = random.randrange(1,10)
# fill upper and lower triangles symmetrically by replicating diagonally
for v in range(1,dim):
iterations = dim - v
x = v
y = 0
while iterations > 0:
new_square[x][y] = new_square[y][x] = random.randrange(1,10)
x += 1
y += 1
iterations -= 1
return new_square
# sanity test
def test_symmetry(square):
dim = len(square[0])
isSymmetric = ''
for x in range(0, dim):
for y in range(0, dim):
if square[x][y] != square[y][x]:
isSymmetric = 'NOT'
print "Matrix is", isSymmetric, "symmetric"
def showSquare(square):
# Print out square matrix
columnHeader = ' '
for i in range(len(square)):
columnHeader += ' ' + str(i)
print columnHeader
i = 0;
for col in square:
print i, col # print row number and data
i += 1
def myMain(argv):
if len(argv) == 1:
nodeCount = 6
else:
try:
nodeCount = int(argv[1])
except:
print "argument must be numeric"
quit()
# keep nodeCount <= 9 to keep the cost matrix pretty
costMatrix = fillCostMatrix(nodeCount)
print "Cost Matrix"
showSquare(costMatrix)
test_symmetry(costMatrix) # sanity test
if __name__ == "__main__":
import sys
myMain(sys.argv)
# vim:tabstop=8:shiftwidth=4:expandtab
Это тривиально для Pythonically заполнения [i][j], если заполняется [j][i]. Вопрос о хранении немного интереснее. Можно увеличить класс массива numpy с помощью атрибута packed, который полезен как для сохранения хранилища, так и для последующего чтения данных.
class Sym(np.ndarray):
# wrapper class for numpy array for symmetric matrices. New attribute can pack matrix to optimize storage.
# Usage:
# If you have a symmetric matrix A as a shape (n,n) numpy ndarray, Sym(A).packed is a shape (n(n+1)/2,) numpy array
# that is a packed version of A. To convert it back, just wrap the flat list in Sym(). Note that Sym(Sym(A).packed)
def __new__(cls, input_array):
obj = np.asarray(input_array).view(cls)
if len(obj.shape) == 1:
l = obj.copy()
p = obj.copy()
m = int((np.sqrt(8 * len(obj) + 1) - 1) / 2)
sqrt_m = np.sqrt(m)
if np.isclose(sqrt_m, np.round(sqrt_m)):
A = np.zeros((m, m))
for i in range(m):
A[i, i:] = l[:(m-i)]
A[i:, i] = l[:(m-i)]
l = l[(m-i):]
obj = np.asarray(A).view(cls)
obj.packed = p
else:
raise ValueError('One dimensional input length must be a triangular number.')
elif len(obj.shape) == 2:
if obj.shape[0] != obj.shape[1]:
raise ValueError('Two dimensional input must be a square matrix.')
packed_out = []
for i in range(obj.shape[0]):
packed_out.append(obj[i, i:])
obj.packed = np.concatenate(packed_out)
else:
raise ValueError('Input array must be 1 or 2 dimensional.')
return obj
def __array_finalize__(self, obj):
if obj is None: return
self.packed = getattr(obj, 'packed', None)
`` `