1
$\begingroup$

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

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – Brian Borchers Sep 26 at 18:23
  • $\begingroup$ There’s an SO post about this: stackoverflow.com/questions/1437501/… $\endgroup$ – Amit Hochman Sep 27 at 7:55
  • $\begingroup$ @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 $\endgroup$ – Bananach Sep 27 at 8:12
  • $\begingroup$ Isn’t math.net, mentioned in the post, adequate? For example, for QR, you have this: collective2.com/c2explorer_help/html/… $\endgroup$ – Amit Hochman Sep 27 at 8:18
  • $\begingroup$ Or this: numerics.mathdotnet.com/api/… $\endgroup$ – Amit Hochman Sep 27 at 8:21
2
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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? $\endgroup$ – Bananach Sep 27 at 12:20
  • $\begingroup$ 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. $\endgroup$ – rchilton1980 Sep 27 at 13:26
  • $\begingroup$ 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 $\endgroup$ – Bananach Oct 4 at 21:28
  • $\begingroup$ 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) $\endgroup$ – rchilton1980 Oct 7 at 0:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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