Search code examples
c++numpyeigen

what is equal to np.linalg.solve(A, B) of Python in eigen C++


#!/usr/bin/python
x = np.linalg.solve(A, B)

The above can calculate the root of Ax = B, here, A is a 3 by 3 matrix and B is a 3 by 1 vector. And I want to find a function in Eigen Library with the same function, instead of Python Numpy Library. The result x is right using Python Numpy linalg.solve(). The matrix A and B is shown below:

A: 
  64  256 1024
  48  256 1280
  24  192 1280
B: 
-9
 0
 0

However, I choose the following code (C++ Eigen) to solve the same problem, the error is shown to me.

// C++ with Eigen Library
auto x = A.colPivHouseholderQr().solve(B)

The above code has an error in running time:

frenet: /usr/include/eigen3/Eigen/src/QR/ColPivHouseholderQR.h:546: void
 Eigen::ColPivHouseholderQR<MatrixType>::_solve_impl(const RhsType&, 
DstType&) const [with RhsType = Eigen::Matrix<double, -1, 1>; DstType = 
Eigen::Matrix<double, -1, 1>; _MatrixType = Eigen::Matrix<double, -1,
 -1>]: Assertion `rhs.rows() == rows()' failed.

I don't know what happens. Wish you could help me soon!


Solution

  • If A is known at compiletime to be 3x3, I'd recommend calculating the inverse directly (only once if you need it multiple times)

    #include <Eigen/LU>
    #include <iostream>
    
    int main() {
      Eigen::Matrix3d A;
      A << 64, 256, 1024, 48, 256, 1280, 24, 192, 1280;
      Eigen::Vector3d B;
      B << -9.0, 0, 0;
      Eigen::Matrix3d A_inv = A.inverse();
      Eigen::Vector3d x = A_inv * B;
    
      std::cout << "solution x=\n" << x << "\n\nresidual A*x-B=\n" << A * x - B << '\n';
    }
    

    For larger A you need to pick one of the decompositions which best suits your problem.

    Also, be very careful with using auto together with Eigen expressions: https://eigen.tuxfamily.org/dox-devel/TopicPitfalls.html#TopicPitfalls_auto_keyword