Numpy Hermitian Matrix class

Question

Are you aware of something like a hermitian matrix class in numpy? I'd like to optimize matrix calculations like

B = U * A * U.H

, where A (and thus, B) are hermitian. Without specification, all matrix elements of B are calculated. In fact, it should be able to save a factor of about 2 here. Do I miss something?

The method I need should take take the upper/lower triangle of A, the full matrix of U and return the upper/lower triangle of B.

Answer 1

I don't think there exists a method for your specific problem, but with a little thought you might be able to build an algorithm from the low-level BLAS routines that are wrapped in SciPy. For example, dgemm, dsymm, and dtrmm do general, symmetric, and triangular matrix products respectively. Here's an example of using them:

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))

There are many other low-level BLAS routines available; I find the NETLIB page to be a good resource to learn what they do. You may be able to cleverly use some combination of the available routines to efficiently solve the problem you have in mind.

Edit: it looks like there are LAPACK routines that quickly compute exactly what you want: dsytrd or zhetrd, but unfortunately these don't appear to be wrapped directly in scipy.linalg.lapack, though scipy does provide cython wrappers for them. Best of luck!

Answer 2

I needed tridiagonal reduction of a symmetric/Hermitian matrix A,

T = Q^H * A * Q

– presumably OP's underlying problem – and I've just submitted a pull request to SciPy for properly interfacing LAPACK's {s,d}sytrd (for real symmetric matrices) and {c,z}hetrd (for Hermitian matrices). All routines use either only the upper or the lower triangular part of the matrix.

Once this has been merged, it can be used like

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

The vectors d and e now contain the main diagonal and the upper and lower diagonal, respectively.