Solving a Least Squares Problem

Here is a simple program which reads in a matrix and a right side, builds a least squares factorization object, and then calculates the solution to the least squares problem. The program is nearly identical to the first factorization example problem in the previous chapter.

#include <iostream> // 1

#include <rw/math/genmat.h>

#include <rw/math/mathvec.h>

#include <rw/lapack/lsqr.h>

#include <rw/lapack/qrcalcp3.h>

int main()

{

RWGenMat<float> A; // 2

RWMathVec<float> b;

std::cin >> A >> b;

RWLeastSqQR<float, RWQRCalcP3<float> > ls(A); // 3

RWMathVec<float> x = ls.solve(b); // 4

std::cout << x << std::endl;

return 0;

}

// 1 The header file lsqr.h includes the declaration of the least squares class RWLeastSqQR<T,QRCalc>. This is the QR decomposition implementation of the least squares class. Other implementations are described subsequently.

// 2 Here the matrix A and the vector b are defined and read from standard input.

// 3 The least squares decomposition is created from the matrix A. The second template parameter for the class RWLeastSqQR indicates that we are using the class RWQRCalcP3<float> to compute the QR decomposition of A.

// 4 The solution vector is calculated. If the system is overdetermined, this will be the solution that best fits in a least squares sense. If the system is underdetermined, this will be the solution with the smallest norm.