2
$\begingroup$

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?

Improve this question
$\endgroup$
8
  • 3
    $\begingroup$ 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 $\endgroup$
    – vibe
    Oct 14, 2020 at 3:54
  • $\begingroup$ 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. $\endgroup$
    – Charlie S
    Oct 14, 2020 at 3:59
  • $\begingroup$ @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. $\endgroup$
    – Vittore Scolari
    Oct 14, 2020 at 9:24
  • 1
    $\begingroup$ @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!) $\endgroup$
    – Charlie S
    Oct 14, 2020 at 13:01
  • 1
    $\begingroup$ 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). $\endgroup$
    – Charlie S
    Oct 14, 2020 at 17:58

1 Answer 1

Reset to default
2
$\begingroup$

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.

Improve this answer
$\endgroup$
2
  • 2
    $\begingroup$ 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 $\endgroup$
    – Vittore Scolari
    Oct 14, 2020 at 9:06
  • 3
    $\begingroup$ Unfortunately, [dsyevd] and friends combine householder diagonalization with givens QR iteration, so they must accumulate them together into tabulated storage. The process you describe (diagonalizing M, transforming A, then diagonalizing it) sounds suspiciously like a generalized eigendecomposition but it's hard to say with certainty. This problem can be especially tidy if either operand is positive definite, like your A, see [dsygv]. If you could post another question, that lays out your larger problem (at the PDE level or semi-discrete level) it might lead to more efficient solutions. $\endgroup$
    – rchilton1980
    Oct 14, 2020 at 13:12

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.