Numpy Hermitian 矩阵类答案

【问题标题】：Numpy Hermitian Matrix classNumpy Hermitian 矩阵类
【发布时间】：2016-02-07 13:57:10
【问题描述】：

你知道像 numpy 中的 Hermitian 矩阵类这样的东西吗？我想优化矩阵计算，如

B = U * A * U.H

，其中 A（因此，B）是厄米特。如果没有指定，则计算 B 的所有矩阵元素。其实这里应该可以节省2倍左右的系数。我错过了什么吗？

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

【问题讨论】：

你这样做有什么问题
好吧，最后我需要三角矩阵 B[numpy.triu_indices(dim),0]。但是，在计算 (UAU.H)[numpy.triu_indices(dim)] 时，我首先计算 (UAU.H) 的所有元素，然后选择上三角形。只计算上三角不是更合理吗？
这里定义的厄米特函数：docs.scipy.org/doc/numpy-dev/user/…
对不起，我想我不够清楚。我不需要计算矩阵的厄米特函数。我想利用自伴随矩阵（或实数的对称矩阵）的对称性，以实现更有效的乘积/转换计算。用更简单的话来说：计算 B 的上三角形而不先计算下三角形。我想避免 python 循环。

标签： performance numpy matrix linear-algebra multiplication

【解决方案1】：

我需要对对称/厄米矩阵 A 进行三对角约简，

T = Q^H * A * Q

- 大概是 OP 的潜在问题 - 我刚刚提交了 a pull request to SciPy 以正确连接 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 现在分别包含主对角线和上下对角线。

【讨论】：

【解决方案2】：

我认为不存在针对您的特定问题的方法，但稍加考虑，您可能能够从包装在 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 例程可用；我发现the NETLIB page 是一个很好的资源来了解他们的工作。您也许可以巧妙地使用一些可用例程的组合来有效地解决您想到的问题。

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

【讨论】：