# Efficient change of basis real positive definite symmetric matrix

I need to optimize a code where the most performance critical part is doing a 'change of basis', in other words it is an unitary similarity transformation on a big real positive definite symmetric matrix real matrix. This consists in the following operation: $$U^T A U$$, with $$A$$ real positive definite symmetric matrix and $$U$$ real unitary.

At the moment, I am achieving this using BLAS DGEMM two times. But I am not very satisfied since this is ignoring that the left hand of $$A$$ is equal to the transpose of the right hand of $$A$$. Also it is ignoring all nice $$A$$ properties.

Looking at all LAPACK routines that do unitary similarity transformations, not a single one seems to actually use DGEMM, am I missing a simple optimization opportunity?

• If $A$ is symmetric, you should use DSYMM instead of DGEMM for the first product. You might also consider a Cholesky factorization of $A$ followed by TRMM and DSYRK
– vibe
Oct 14 '20 at 3:54
• Is your matrix sparse? Intel MKL has a routine called mkl_sparse_sypr (or mkl_sparse_?_syprd if you want a dense result) that performs the product $U^T A U$, where $A$ must be a symmetric matrix, and returns the upper triangular portion of the result. This is way faster than any combination of routines that I have tested. Unfortunately, I haven't seen this in any other linear algebra libraries. Oct 14 '20 at 3:59
• @vibe I am considering the first option, it may gain a fraction of speed but the algorithm and data structures stays the same. Regarding the Cholesky factorization, can you please add some details? Maybe even a full answer? I am interested since the code is doing Cholesky in another place for another purpose and I might save it. Oct 14 '20 at 9:24
• @VittoreScolari mkl_sparse_?_syprd takes a dense $U$. Your $A$ is essentially ~50% sparse, so I think its worth a shot. I have tried @vibe 's approach as well: DSYRK performs the operation $B^T B$, where in your case $B = L^T U$. I found syprd to be faster (your experience may be different!) Oct 14 '20 at 13:01
• A similar approach without cholesky: $A = L + L^T$ (an additive decomposition). $B = LU$, $U^T A U = (U^T B) + (U^T B)^T$. It involves a trivial decomposition, one TRMM, one DGEMM, and some addition (you only need to compute your preferred triangle, as the result is symmetric). Not sure if this is better in practice as I have never tried it. I have a feeling Intel does something like this because they don't require $A$ to be positive definite (hence no internal cholesky). Oct 14 '20 at 17:58

For the most part, LAPACK represents unitary matrices like $$\mathbf U$$ as products of householder reflectors, and provides specialized routines to work with that representation (for instance, use [dormqr] to multiply by $$\mathbf U$$, or [dorgqr] to explicitly tabulate $$\mathbf U$$ into a dense matrix).
If you already have $$\mathbf U$$ in householder format, it would be more idiomatic to use [dormqr] to update A, instead of [dgemm]. In particular, two calls to [dormqr] can update $$\mathbf B = \mathbf U^T \mathbf A \mathbf U$$ in place (overwriting $$\mathbf A$$ with $$\mathbf B$$), while [dgemm] will require a temporary matrix. This point might be somewhat academic, though, as [dormqr] requires additional workspace while [dgemm] does not. Although there is a workspace-free routine [dorm2r] that performs a similar function as [dormqr], it's not recommended because it's a BLAS2 algorithm and won't be as fast.
I suppose if I already was just given $$\mathbf U$$ explicitly tabulated, I'd probably just stick with [dgemm]. But if you have control over how $$\mathbf U$$ is generated/computed (for instance, some hand-rolled gram-schmidt procedure applied to some other basis set), I'd strongly consider refactoring everything to use the LAPACK tooling instead (that is use [dgeqrf] to orthogonalize your basis into $$\mathbf U$$ and then use [dormqr] to apply it, there might not be any need to ever form it explicitly).
Out of curiousity, what do you do with $$\mathbf B$$ once you have it? There might be more refactoring/optimization opportunities if you can share more about the larger process.
• The way I get $U$ is by diagonalising another matrix $M$ through [dsyevd]. I wonder if there is a way to get reflectors instead of vectors from the diagonalization procedure, and if this can help in performances? Regarding your second question, I use the $B$ diagonal to solve a PDE problem, but I really doubt I can optimize that part of the code since it is fast Oct 14 '20 at 9:06