Hello, I have weighted least squares problem. This problem can be solved by famous equation:
W = (XT * Z * X)^-1 * (XT * Z * Y)
and simpler version W = (XT * X)^-1 * (XT *Y)
where
XT = X.transpose()
W = ((XT * Z * X).inverse()) * (XT * Z * Y)
Sizes (rows, cols) of matrices are:
X - (m, n)
Z - (m, m)
Y - (m, 1)

This equation is in core of my system, so, I need it be calculated as fast as possible. I'm beginner in Eigen, so I need your advice - few lines of C++ code, which will solve weighted least squares problem as fast as possible using Eigen.

Conditions:
1. My IDE is Visual Studio 2012.
2. My CPU supports AVX (Core I7 Sandy Bridge), but I also need executable which relies on SSE4.2 only (without AVX).
3. My matrix sizes are m < 1000, n < 100.
4. Matrix Z is diagonal.

Your code should be connected with my code:

double *matrixX = new double [m*n];
double *matrixZ = new double [m*m]; // or, I can define double *matrixZ = new double [m*1];
double *matrixY = new double [m];
// Here I defined my matrices.

// Here I use your C++ code

Questions
1. How should I set up environment to achieve maximum speed?
2. What I have to do if matrix (XT * Z * X) is singular? How to determine singularity?
3. I need your C++ code with Eigen.

Hi garoldq,

For the first part of your question, there is general information here: https://eigen.tuxfamily.org/dox-devel/TopicMultiThreading.html
Also, if you can (for example if you have to solve your least square problem many times), you can parallelize your code (eg with OpenMP).

There is the function Map (documented here: https://eigen.tuxfamily.org/dox/classEigen_1_1Map.html) to construct vectors from vectors of double.
As an example, the following code copy a vector<double> into the first column of a matrix.

Code: Select all

MatrixXd m;
m.resize(3,2);
   
vector< double > v( { 1.1, 1.2, 1.3} );
m.col(0) = VectorXd::Map( &v[0], 0, 3 );


Concerning your question 2., I suppose the answer depends on your problem and not on Eigen.

Probably, members of this forum won't write you a complete piece of code but I'm happy to help you if I can.

Best,

Raf

Some ressources for LS solving in Eigen: http://eigen.tuxfamily.org/dox-devel/gr ... uares.html

and a benchmark of the different choices to help you pick the best tradeoff:

http://eigen.tuxfamily.org/dox-devel/gr ... hmark.html

if perf matters the most, start with the normal-equation + LLT (or LDLT if small enough), check the status if the dec by looking at the diagonal entries, if the ratio between largest and smallest entry is too large, then fall back to CompleteOrthogonalDecomposition for the minimal norm solution.