I would like to solve the linear system Ax = b
in a linear least squares fashion, thereby obtaining x
. Matrices A
, x
and b
contain elements that are complex numbers.
Matrix A
has dimensions of n
by n
, and A
is a square matrix that is also lower triangular. Vectors b
and x
have lengths of n
. There are as many unknowns as there are equations in this system, but since b
is a vector filled with actual measured "data", I suspect that it would be better to do this in a linear least squares fashion.
I am looking for an algorithm that will efficiently solve this system in a LLS fashion, using perhaps a sparse matrix data structure for lower-triangular matrix A
.
Perhaps there is a C/C++ library with such an algorithm already available? (I suspect that it is best to use a library due to optimized code.) Looking around in the Eigen matrix library, it appears that SVD decomposition can be used to solve a system of equations in a LLS fashion (link to Eigen documentation). However, how do I work with complex numbers in Eigen?
It appears that the Eigen library works with the SVD, and then uses this for LLS solving.
Here is a code snippet demonstrating what I would like to do:
#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
Here is the compiler error:
error: missing template arguments before '(' token
UPDATE
Here is an updated program showing how to deal with the LLS solving using Eigen. This code does indeed compile correctly.
#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