Using LAPACK to compute $B^{-1}AB^{-T}$ for thin $B$

Question

How can I use BLAS/LAPACK to compute $$ B^{-1}AB^{-T} $$ where $A\in\mathbb{R}^{n,n}$, $B\in\mathbb{R}^{m,n}$ is full rank matrix with $m>n$, and $B^{-1}y:=\arg \min_{x} \|Bx-y\|_{2}$.

In theory, one could compute $$ B=QR $$ with orthogonal $Q\in\mathbb{R}^{m,n}$ and triagonal $R\in\mathbb{R}^{n,n}$ and then $$ R^{-1}Q^{T}AQR^{-T}. $$

However, I don't know how to implement this in BLAS/LAPACK (in C# if it matters). I haven't used either and the documentations don't help me at all. I found geqrf but have no idea what I should do with the outputs it gives me.

Alternatively: does anyone know of a widely available higher-level library that can do this for me in C#?

Answer 1

You should be able to do this efficiently using LAPACK/BLAS, using the QR factorization [geqrf], orthogonal multiplication [ormqr], and triangular solve [trsm] routines. LAPACK and BLAS place a lot of emphasis on economy of storage, so many of these algorithms operate in-place, overwriting their inputs with their outputs (geqrf, ormqr and trsm all work like this).

You should start off with geqrf, it will overwrite $\mathbf B$ with both $\mathbf Q$ and $\mathbf R$. The $\mathbf R$ output is just tabulated directly into the upper triangle, but $\mathbf Q$ is represented compactly in "reflector form". This representation requires a small/vector-sized amount of additional workspace ($\mathbf \tau$), empty storage that you provide to geqrf, that it will fill upon exit. You need to hang onto $\mathbf \tau$ to do anything else with $\mathbf Q$.

Because $\mathbf Q$ is stored in that special/reflector form, you can't just use BLAS/gemm to apply it. Instead you should use the LAPACK routine ormqr, it can apply $\mathbf Q$ to $\mathbf A$ from either the left or right and also apply transposition. It also operates inplace, overwriting $\mathbf A$ with (eg) $\mathbf Q^T \mathbf A$ or $\mathbf A \mathbf Q$. This can make it a little tricky to reason about, because (eg) $\mathbf Q^T \mathbf A$ and $\mathbf A$ are not necessarily the same size ($\mathbf Q$ being tall/skinny). Nevertheless you should be able to use two calls to ormqr (from the 'L'eft with a 'T'ranspose, then from the 'R'ight with 'N'otranspose) to get $\mathbf C = \mathbf Q^T\mathbf A\mathbf Q$ overwritten into the upper left corner of your original $\mathbf A$.

The last step, forming $\mathbf R ^{-1} \mathbf C \mathbf R^{-T}$ can be accomplished using the BLAS routine trsm. It applies the inverse of a triangular system, overwriting the input $\mathbf C$ with (eg) the output $\mathbf R^{-1} \mathbf C$. It can apply $\mathbf R$ from the left or right, and also apply a transposition to it. You'll need two calls, one from the 'L'eft with 'N'otranspose, then another from the 'R'ight with 'T'ranspose. In both cases you'll want to use uplo='U'pper because that's where geqrf stores $\mathbf R$.

Unfortunately I am not a C# expert, but I would imagine there's some usual way to call native/unmanaged "extern C" libraries, that's where you'd want to start.