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

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#?

• This would be relatively easy to implement in LAPACK, but "How do I do this in LAPACK?" questions by themselves are off-topic for this group. You could formulate this question by backing up and explaining the larger problem that you're working on. – Brian Borchers Sep 26 '19 at 18:23
• @Amit I already have a Blast/Lapack library, but I don't know how to use it. I am used to numpy/scipy levels of abstraction – Bananach Sep 27 '19 at 8:12
• Isn’t math.net, mentioned in the post, adequate? For example, for QR, you have this: collective2.com/c2explorer_help/html/… – Amit Hochman Sep 27 '19 at 8:18
• Or this: numerics.mathdotnet.com/api/… – Amit Hochman Sep 27 '19 at 8:21

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.

• Works perfectly. Just one more question. To accelerate the second call to ormqr I guess you should create a new array that contains only the first few matrix rows of the overwritten A (i.e. a new array that contains only Q^TA, nothing else) so that the call takes . Is there an elegant way to do this? – Bananach Sep 27 '19 at 12:20
• Yes, you can/should do this. You don't even need to make a copy, you can just "slice" top(A) by adjusting the inputs to ormqr that specify the size/layout of A. [Specifically, top(A) has the same number of columns, leading dimension and data pointer as A, it just has fewer rows]. When implemented correctly, you should see that bottom(A) is untouched by the second call to ormqr. – rchilton1980 Sep 27 '19 at 13:26
• Are you sure this is also possible in a column major BLAS implementation? The top matrix isn't a contiguous block of memory in that case – Bananach Oct 4 '19 at 21:28
• Yes, should be OK. You'll just need to adjust the rows/columns arguments but leave the leading dimension alone (since you have LD != rows, LAPACK will "jump" across the bottom parts of A) – rchilton1980 Oct 7 '19 at 0:08