CODEWITHSUNDEEP
c++numpyeigen
Bryan YU
Bryan YU

Reputation: 77

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!

Upvotes: 0

Views: 2673

Answers (1)

chtz
chtz

Reputation: 18827

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

Upvotes: 1

Related Questions