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$
  • $\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$ – Christian Clason Sep 17 '16 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 Sep 17 '16 at 13:07
  • $\begingroup$ Yes, real or complex shouldn't make a difference here. $\endgroup$ – Christian Clason Sep 17 '16 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 Sep 17 '16 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$ – Christian Clason Sep 17 '16 at 14:06
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$
  • $\begingroup$ The internal implementation for A*B.inverse() is really inverse B first and then do the gemm? $\endgroup$ – lorniper Sep 17 '16 at 15:14
  • $\begingroup$ Yes, in C++ the precedence of the dot(.) operator is much higher than *. $\endgroup$ – Bill Greene Sep 17 '16 at 15:31

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.