Numpy Hermitian矩陣類

Question

你知道像numpy中的厄密矩陣類的東西嗎？我想，以優化矩陣計算等

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程序可用;我發現the NETLIB page是一個很好的資源，可以學習他們的工作。您可以巧妙地使用可用例程的一些組合來有效地解決您所考慮的問題。

編輯：它看起來像有LAPACK例程，快速的計算正是你想要的：dsytrd或zhetrd，但不幸的是，這些似乎並沒有被在scipy.linalg.lapack直接包裹，儘管SciPy的確實提供cython wrappers他們。祝你好運！

Answer 2

我需要三對角縮小對稱的/埃爾米特矩陣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現在包含在主對角線和上部和下部對角，分別。