Ax = b
線形最小二乗法で線形システムを解き、それによってを取得したいと思いx
ます。行列A
、x
およびb
複素数である要素が含まれています。
行列A
の次元はn
byn
であり、A
これも下三角行列です。ベクトルb
とx
の長さはn
。このシステムには方程式と同じ数の未知数がありますが、b
は実際に測定された「データ」で満たされたベクトルであるため、線形最小二乗法でこれを行う方がよいと思います。
私は、おそらく下三角行列のスパース行列データ構造を使用して、LLS方式でこのシステムを効率的に解決するアルゴリズムを探していますA
。
おそらく、そのようなアルゴリズムがすでに利用可能なC / C ++ライブラリがありますか?(コードが最適化されているため、ライブラリを使用するのが最適だと思います。)Eigen行列ライブラリを見ると、SVD分解を使用してLLS方式で連立方程式を解くことができるようです(Eigenドキュメントへのリンク)。しかし、Eigenで複素数を操作するにはどうすればよいですか?
EigenライブラリはSVDで動作し、これをLLS解法に使用しているようです。
これが私がやりたいことを示すコードスニペットです:
#include <iostream>
#include <Eigen/Dense>
#include <complex>
using namespace Eigen;
int main()
{
// I would like to assign complex numbers
// to A and b
/*
MatrixXcd A(4, 4);
A(0,0) = std::complex(3,5); // Compiler error occurs here
A(1,0) = std::complex(4,4);
A(1,1) = std::complex(5,3);
A(2,0) = std::complex(2,2);
A(2,1) = std::complex(3,3);
A(2,2) = std::complex(4,4);
A(3,0) = std::complex(5,3);
A(3,1) = std::complex(2,4);
A(3,2) = std::complex(4,3);
A(3,3) = std::complex(2,4);
*/
// The following code is taken from:
// http://eigen.tuxfamily.org/dox/TutorialLinearAlgebra.html#TutorialLinAlgLeastsquares
// This is what I want to do, but with complex numbers
// and with A as lower triangular
MatrixXf A = MatrixXf::Random(3, 3);
std::cout << "Here is the matrix A:\n" << A << std::endl;
VectorXf b = VectorXf::Random(3);
std::cout << "Here is the right hand side b:\n" << b << std::endl;
std::cout << "The least-squares solution is:\n"
<< A.jacobiSvd(ComputeThinU | ComputeThinV).solve(b) << std::endl;
}// end
コンパイラエラーは次のとおりです。
error: missing template arguments before '(' token
アップデート
これは、Eigenを使用してLLS解決に対処する方法を示す更新されたプログラムです。このコードは確かに正しくコンパイルされます。
#include <iostream>
#include <Eigen/Dense>
#include <complex>
using namespace Eigen;
int main()
{
MatrixXcd A(4, 4);
A(0,0) = std::complex<double>(3,5);
A(1,0) = std::complex<double>(4,4);
A(1,1) = std::complex<double>(5,3);
A(2,0) = std::complex<double>(2,2);
A(2,1) = std::complex<double>(3,3);
A(2,2) = std::complex<double>(4,4);
A(3,0) = std::complex<double>(5,3);
A(3,1) = std::complex<double>(2,4);
A(3,2) = std::complex<double>(4,3);
A(3,3) = std::complex<double>(2,4);
VectorXcd b(4);
b(0) = std::complex<double>(3,5);
b(1) = std::complex<double>(2,0);
b(2) = std::complex<double>(8,2);
b(3) = std::complex<double>(4,8);
std::cout << "Here is the A matrix:" << std::endl;
std::cout << A << std::endl;
std::cout << "Here is the b vector:" << std::endl;
std::cout << b << std::endl;
std::cout << "The least-squares solution is:\n"
<< A.jacobiSvd(ComputeThinU | ComputeThinV).solve(b) << std::endl;
}// end