2

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.

ali_m
  • 71,714
  • 23
  • 223
  • 298
newbie
  • 51
  • 6
  • what problem are you facing doing this – Kumar Saurabh Nov 06 '15 at 12:09
  • Well, in the end I need the triangular matrix B[numpy.triu_indices(dim),0]. However, calculating (U*A*U.H)[numpy.triu_indices(dim)], I do first calculate all elements of (U*A*U.H) and then pick the upper triangle. Wouldnt it be more reasonable to calculate the upper triangle only? – newbie Nov 06 '15 at 12:32
  • Hermitian function defined here: https://docs.scipy.org/doc/numpy-dev/user/numpy-for-matlab-users.html#customizing-your-environment – Lee Nov 06 '15 at 13:13
  • Sorry, I think i was not clear enough. I do not need a function that calculates the hermitian of a matrix. I'd like to exploit the symmetry of self-adjoint matrices (or symmetric matrices for real numbers) in order to achive a more efficient calculation of products / transformations. In more easy words: Calculate the upper triangle of B without calculating the lower triangle first. I'd like to avoid python loops. – newbie Nov 06 '15 at 13:25

2 Answers2

4

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!

jakevdp
  • 77,104
  • 11
  • 125
  • 160
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.

Nico Schlömer
  • 53,797
  • 27
  • 201
  • 249