Deb*_*ski 66 python numpy matrix
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
一个自动的Hermitian也会很好,虽然在写作的时候我不需要它.
Eri*_*got 70
如果您能够在计算之前对矩阵进行对称化,则以下内容应该相当快:
def symmetrize(a):
return a + a.T - numpy.diag(a.diagonal())
这在合理的假设下工作(例如在运行之前不做这两者a[0, 1] = 42和矛盾).a[1, 0] = 123symmetrize
如果你真的需要一个透明的对称化,你可以考虑继承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之前必须稍微调整代码:def __setitem__(self, indexes, value): (i, j) = indexes...
Mat*_*att 20
numpy中对称矩阵的最优处理的更一般问题也让我感到烦恼.
在研究之后,我认为答案可能是numpy受到对称矩阵的底层BLAS例程支持的内存布局的限制.
虽然一些BLAS例程确实利用对称性来加速对称矩阵的计算,但它们仍然使用与完整矩阵相同的存储器结构,即n^2空间而不是n(n+1)/2.只是他们被告知矩阵是对称的,并且只使用上三角或下三角中的值.
一些scipy.linalg例程接受传递给BLAS例程的标志(如sym_pos=Trueon linalg.solve),虽然numpy中对此的更多支持会很好,特别是像DSYRK(对称等级k更新)这样的例程的包装器,这将允许Gram矩阵计算比dot(MT,M)快一点.
(可能会因为担心在时间和/或空间上优化2倍常数因素而感到挑剔,但它可以对您在单个机器上管理的问题的大小有所影响......)
小智 7
存在许多众所周知的存储对称矩阵的方法,因此它们不需要占用n ^ 2个存储元件.此外,重写共同操作以访问这些修订的存储方式是可行的.最权威的工作是Golub和Van Loan,Matrix Computations,1996年第3版,约翰霍普金斯大学出版社,第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索引.
Golub和Van Loan提供了一个算法1.2.3,它展示了如何访问这样存储的V来计算y = V x + y.
Golub和Van Loan还提供了一种以对角线显性形式存储矩阵的方法.这不会节省存储空间,但支持对某些其他类型的操作进行随时访问.