QR_lssolve() seems to insist on only solving square problems, unlike gsl_linalg_QR_lssolve()
crazzell opened this issue · 2 comments
I am trying to use gsl (1.16.0.6) for solving over-determined least-squares problems, and a simple test case is not working for me.
require("gsl")
include GSL
m = GSL::Matrix::alloc([0.18, 0.60, 0.57], [0.24, 0.99, 0.58],
[0.14, 0.30, 0.97], [0.51, 0.19, 0.85], [0.34, 0.91, 0.18])
B = GSL::Vector[1, 2, 3, 4, 5]
qr, tau = m.QR_decomp
x, r = qr.QR_lssolve(tau,B)
Whereas, if I write the same problem in plain C, it works find:
#include <stdio.h>
#include <gsl/gsl_linalg.h>
int
main (void)
{
double a_data[] = {0.18, 0.60, 0.57, 0.24, 0.99, 0.58,
0.14, 0.30, 0.97, 0.51, 0.19, 0.85, 0.34, 0.91, 0.18};
double b_data[] = {1.0,2.0,3.0,4.0,5.0};
gsl_matrix_view m = gsl_matrix_view_array (a_data, 5, 3);
gsl_vector_view b = gsl_vector_view_array (b_data, 5);
gsl_vector *x = gsl_vector_alloc (3); // size equal to n
gsl_vector *residual = gsl_vector_alloc (5); // size equal to m
gsl_vector *tau = gsl_vector_alloc (3); //size equal to min(m,n)
gsl_linalg_QR_decomp (&m.matrix, tau); //
gsl_linalg_QR_lssolve(&m.matrix, tau, &b.vector, x, residual);
printf ("x = \n");
gsl_vector_fprintf (stdout, x, "%g");
gsl_vector_free (x);
gsl_vector_free (tau);
gsl_vector_free (residual);
return 0;
}
x =
8.06835
0.884443
0.231883
The only way I can get the ruby code to work is by making the problem square, but that defeats the overall purpose of a least-squares fit in my use-case. Of course, it is just possible that I am messing up the row and column dimensions in the Ruby code in a way that I can't see for looking, but I simply can't find another permutation that should work. I am guessing the ruby-wrapper is mistakenly assuming that over-determined problems are always user-error, whereas this is precisely where a LS fit is most useful.
Thanks! Charles.
Forgot to add the error message from Ruby with the above script:
../testls.rb:9:in 'QR_lssolve': Ruby/GSL error code 19, matrix size must match solution size (file qr.c, line >193), matrix/vector sizes are not conformant (GSL::ERROR::EBADLEN)
'
from ../testls.rb:9:in '
Apologies for the false alarm, I found that the following method works:
require("gsl")
include GSL
m = GSL::Matrix::alloc([0.18, 0.60, 0.57], [0.24, 0.99, 0.58],
[0.14, 0.30, 0.97], [0.51, 0.19, 0.85], [0.34, 0.91, 0.18])
B = GSL::Vector[1, 2, 3, 4, 5]
qr, tau = m.QR_decomp
x = GSL::Vector.alloc(3)
r = GSL::Vector.alloc(5)
qr.QR_lssolve(tau,B,x,r)
p x
I was not expecting to have to put the outputs {x,r} as arguments to the function, so this took me a while to find.