Is there an Eigen C ++ library equivalent to gsl_multifit_wlinear ()?
I'm trying to port some old GSL code to use only the Eigen header. I do this as a side job for some data scientists, and so my math has gone a bit out of college. I am trying to see if there is an equivalent function or module, or even a more general linear algebra term for calling the gsl_multifit_wlinear () function. It is similar to the least squares curve fitting function.
If there is a way to directly convert, what would it be?
If not, is there any other library I can use? Keep in mind that it does not need to be licensed under the GPL or similar "Shared All" license. MIT or BSD is preferred, LGPL and Mozilla / Apache are fine too.
Thank.
source to share
It looks like it gsl_multifit_linear(X,y,c)
solves the problem min_c ||Xc-y||^2
. I'm not really 100% sure from the documentation, but it looks like it gsl_multifit_wlinear(X,y,w,c)
decides min_c ||Xc - y||^2_W
where W = diag(w)
and ||e||^2_W = e'*W^(-1)*e
.
So, you can solve this in Eigen by rewriting min_c ||Xc - y||^2_W
as min_c ||W^(-1/2) (Xc - y)||^2
.
We should:
- compute W ^ (- 1/2), which is simply diagonal with elements W (i, i) = 1 / sqrt (w (i)).
- Calculate A = W ^ (- 1/2) * X and b = W ^ (- 1/2) * y
- Decide
min ||Ac - b||^2
what is the sameAc = b
as the least squares solution
This should work then if you already have Eigen :: Matrix <...> of X, y, w:
Eigen::VectorXd inv_sqrtw = 1.0 / w.array().sqrt();
Eigen::MatrixXd W12 = inv_sqrtw.asDiagonal();
Eigen::MatrixXd A = W12 * X;
Eigen::VectorXd b = W12 * y;
// now solve system
Eigen::VectorXd c = A.jacobiSvd(Eigen::ComputeThinU | Eigen::ComputeThinV).solve(b);
and your answer is in c
More information on solving the least squares problem in Eigen is at http://eigen.tuxfamily.org/dox-devel/group__LeastSquares.html
source to share