Numpy Hermitian Matrix class

Question

您知道numpy中的hermitian矩阵类吗？我想优化矩阵计算，如

B=U*A*U.H

，其中A（因此，B）是厄密的。在没有指定的情况下，计算B的所有矩阵元素。事实上，这里应该能够节省大约2的因子。我错过什么了吗？

我需要的方法应该取A的上/下三角形，U的全矩阵，并返回B的上/下部三角形。

Answer 1

我不认为有一种方法可以解决您的特定问题，但只要稍微考虑一下，您就可以从SciPy中封装的低级BLAS例程构建算法。例如，dgemm、dsymm和dtrmm分别做一般矩阵乘积、对称矩阵乘积和三角矩阵乘积。这里有一个使用它们的例子：

from scipy.linalg.blas import dgemm, dsymm, dtrmm
A = np.random.rand(10, 10)
B = np.random.rand(10, 10)
S = np.dot(A, A.T)  # symmetric matrix
T = np.triu(S)  # upper triangular matrix
# normal matrix-matrix product
assert np.allclose(dgemm(1, A, B), np.dot(A, B))
# symmetric mat-mat product using only upper-triangle
assert np.allclose(dsymm(1, T, B), np.dot(S, B))
# upper-triangular mat-mat product
assert np.allclose(dtrmm(1, T, B), np.dot(T, B))

还有许多其他低级别的BLAS例程可用；我发现NETLIB页面是学习他们做什么的好资源。你可以巧妙地使用一些可用例程的组合来有效地解决你心中的问题。

编辑：看起来有LAPACK例程可以快速计算出您想要的东西：dsytrd或zhetrd，但不幸的是，这些例程似乎没有直接封装在scipy.linalg.lapack中，尽管scipy确实为它们提供了cython封装器。祝你好运！

Answer 2

我需要对称/埃尔米特矩阵A的三对角归约，

T = Q^H * A * Q

——大概是OP的根本问题——我刚刚向SciPy提交了一个pull请求，要求正确连接LAPACK的{s,d}sytrd（对于实对称矩阵）和{c,z}hetrd（对于Hermitian矩阵）。所有例程都只使用矩阵的上三角部分或下三角部分。

一旦合并，就可以像一样使用

import numpy as np
n = 3
A = np.zeros((n, n), dtype=dtype)
A[np.triu_indices_from(A)] = np.arange(1, 2*n+1, dtype=dtype)
# query lwork -- optional
lwork, info = sytrd_lwork(n)
assert info == 0
data, d, e, tau, info = sytrd(A, lwork=lwork)
assert info == 0

向量d和e现在分别包含主对角线以及上对角线和下对角线。