0

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
8
  • 1
    and what line does that error refer to? Dec 6, 2012 at 15:53
  • @RodyOldenhuis: The line refers to A(0,0) = std::complex(3,5);, and the assignments of object type std::complex. Dec 6, 2012 at 15:54
  • Can you show what MatrixXcd is?
    – SinisterMJ
    Dec 6, 2012 at 15:55
  • 1
    Wait, your matrix is already a lower-triangular matrix? Why don't you simply use forward and back substitution?
    – Zeta
    Dec 6, 2012 at 15:57
  • 2
    Wouldn't you have to use std::complex<double>(2,4); instead?
    – SinisterMJ
    Dec 6, 2012 at 16:03

1 Answer 1

3

Since std::complex is a template class, and you init with std::complex(1,1); the compiler doesn't know what type it is.

Use std::complex<double>(1, 1); instead.

1
  • Thank you; you are right about the compiler not knowing the type of the object. Dec 6, 2012 at 16:40

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.