1
$\begingroup$

I know that, from numerical point of view, computing

Ax = b
B=inv(A), x= B*b

are completely different things, and we should factor the matrix using TRF routine then solve with TRS/TRI routine, or in a combined gesv call.

Anyway is there an efficiently routine to specifically deal with

C = A* inv(B)

BTW, I know some templating libraries(e.g. Eigen) might be able to do it well, but I am more interested in a Lapack/MKL solution, that said, correct solution using Eigen is also welcome.

Edit, as Christian Clason pointed out: C = A* inv(B) can be computed via gesv(trans(B),trans(A)), However, in a real implementation with lapack, gesv will overwrite its input, so does it mean we have to introduce temp variable in order to have C = A * inv(B) correctly computed in semantic sense.

$\endgroup$
7
  • $\begingroup$ Note that the columns of $A^{-1}B$ are given by $A^{-1}$ applied to the columns of $B$, so to obtain $X:=A^{-1}B$ you can apply any of the standard methods solve the system $AX=B$ (the LAPACK routines take multiple right hand sides as a matrix, so should work out of the box). For $AB^{-1}$, note that $(AB^{-1})^T = (B^T)^{-1}A^T$. $\endgroup$ Commented Sep 17, 2016 at 10:49
  • $\begingroup$ @ChristianClason I think you made a good point, so we can still use the same lapack routine with complex conjugate in case of complex value? $\endgroup$
    – lorniper
    Commented Sep 17, 2016 at 13:07
  • $\begingroup$ Yes, real or complex shouldn't make a difference here. $\endgroup$ Commented Sep 17, 2016 at 13:38
  • $\begingroup$ @ChristianClason maybe you can leave your answer below, on one hand, make sure "real or complex shouldn't make a difference" means trans() for real and conj() for complex, on the other hand, I am still struggling for a semantically reasonable solution because gesv routines overwrite the input A and B. $\endgroup$
    – lorniper
    Commented Sep 17, 2016 at 13:45
  • $\begingroup$ Maybe you could first edit your question to make it a bit more explicit what you actually want to compute; I was guessing that you really needed the matrix $C=AB^{-1}$ (where $B$ is arbitrary and not, as in your first block, the inverse of $A$) for some reason. If all you need it for is to apply it to some vector(s) $x$ later, it is usually far cheaper to first solve $By=x$ and then compute $Ay$. (Similarly, if you want to solve $Cx=b$, you can do this in two stages, where the first is really cheap since you already have the inverse of $B^{-1}$.) $\endgroup$ Commented Sep 17, 2016 at 14:06

1 Answer 1

1
$\begingroup$

Here is the Eigen example you requested. It follows the approach in Christian Clason's comment and also does the same computation using an explicit inverse for comparison.

#include <iostream>
using std::cout;
using std::endl;
#include <ctime>

#include <Eigen/Core>
#include <Eigen/LU>

void compareInvAndFactor()
{

  typedef Eigen::MatrixXd Matrix;
  const int n = 1000, m = 600;
  Matrix B = Matrix::Random(n, n);
  Matrix A = Matrix::Random(m, n);

  std::clock_t start = std::clock();
  auto Clu = B.transpose().lu().solve(A.transpose()).transpose();
  cout << "Elapsed time using LU factorization = " << (std::clock() - start) / 
          (double)CLOCKS_PER_SEC << " seconds." << endl;

  start = std::clock();
  auto Cinv = A*B.inverse();
  cout << "Elapsed time using inverse = " << (std::clock() - start) / 
          (double)CLOCKS_PER_SEC << " seconds." << endl;

  double diff = (Cinv - Clu).cwiseAbs().maxCoeff();
  cout << "max difference in C matrices =" << diff << endl;
}
$\endgroup$
2
  • $\begingroup$ The internal implementation for A*B.inverse() is really inverse B first and then do the gemm? $\endgroup$
    – lorniper
    Commented Sep 17, 2016 at 15:14
  • $\begingroup$ Yes, in C++ the precedence of the dot(.) operator is much higher than *. $\endgroup$ Commented Sep 17, 2016 at 15:31

Your Answer

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

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