LinAlg SV Doc?

["Kevin H. Hobbs" <hobbsk@ohiou.edu>]

[-- Attachment #1.1: Type: text/plain, Size: 1700 bytes --]

I'm fairly new to using the linear algebra routines so this may be
obvious to everybody else but :

In gsl_linalg_SV_solve how long should b and x be?

Could an example be added for this section?

I believe the length of b must be M and the length of x must be N but I
get really confused here, mostly around :

is M the number of samples or the number of parameters (below)?

I will gladly contribute the function I am working on now as an example.

My problem of the moment is to find the best fitting plane for sets of 7
dimensional points, though 3-D will be enough to get me started.

I have a set of 9 points from f(x,y)=z, and I want to find the best
fitting plane, ax + by + c = z, though these points.

My points are :

  f(1, 1) = 2
  f(1, 2) = 2.11
  f(1, 3) = 3.21
  f(2, 1) = 3.97
  f(2, 2) = 3.76
  f(2, 3) = 4.01
  f(3, 1) = 3.40
  f(3, 2) = 5.47
  f(3, 3) = 4.35

When I put these into the equation for a plane this becomes :

  1a + 1b + 1c = 2
  1a + 2b + 1c = 2.11
  1a + 3b + 1c = 3.21
  2a + 1b + 1c = 3.97
  2a + 2b + 1c = 3.76
  2a + 3b + 1c = 4.01
  3a + 1b + 1c = 3.40
  3a + 2b + 1c = 5.47
  3a + 3b + 1c = 4.35

Then I create the matrix A and the vectors x and b for the linear system
Ax = b (I know x and b are getting confused):

      1 1 1      2.00
      1 2 1      2.11
      1 3 1      3.21
      2 1 1      3.97
  A = 2 2 1  b = 3.76  x = a b c
      2 3 1      4.01
      3 1 1      3.40
      3 2 1      5.47
      3 3 1      4.35


The output I get from the attached program is :

0.983333 0.366667 0.886667

which puts a nice plane through my points.

Do I have my Ms and Ns in the right places?

[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #1.2: find_interpolation_parameters.h --]
[-- Type: text/x-chdr; name="find_interpolation_parameters.h", Size: 455 bytes --]

#ifndef FIND_INTERPOLATION_PARAMETERS_H
#define FIND_INTERPOLATION_PARAMETERS_H
#include <stdio.h>
#include <stdlib.h>
#include <gsl/gsl_linalg.h>

int find_interpolation_parameters( const size_t dim, 
  double * point_data, 
  double * value_data, 
  double * parameters );


int find_interpolation_parameters_SVD( const size_t dim, 
  const size_t num,
  double * point_data, 
  double * value_data, 
  double * parameters );
#endif


[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #1.3: find_interpolation_parameters_SVD.c --]
[-- Type: text/x-csrc; name="find_interpolation_parameters_SVD.c", Size: 1217 bytes --]

#include "find_interpolation_parameters.h"

int find_interpolation_parameters_SVD( const size_t dim, 
  const size_t num,
  double * point_data, 
  double * value_data,
  double * parameters )
{
  gsl_matrix * AU   = gsl_matrix_alloc( num,     dim + 1 );
  gsl_matrix * V    = gsl_matrix_alloc( dim + 1, dim + 1 );
  gsl_vector * S    = gsl_vector_alloc( dim + 1 );
  gsl_vector * work = gsl_vector_alloc( dim + 1 );
  gsl_vector * b    = gsl_vector_alloc( num );
  gsl_vector * x    = gsl_vector_alloc( dim + 1 );

  size_t i;
  size_t j;

  for ( i = 0; i < num; ++i )
    {
    for ( j = 0; j < dim; ++j )
      {
      gsl_matrix_set( AU, i, j, point_data[ i * dim + j ] );
      }
    gsl_matrix_set( AU, i, dim, 1.0 );
    }

  for ( i = 0; i < num; ++i )
    {
    gsl_vector_set( b, i, value_data[i] );
    }

  gsl_linalg_SV_decomp( AU, V, S, work );
  gsl_linalg_SV_solve(  AU, V, S, b, x );

  for ( i = 0; i < dim + 1; ++i )
    {
    parameters[i] = gsl_vector_get( x, i );
    }

  gsl_matrix_free( AU );
  gsl_matrix_free( V );
  gsl_vector_free( S );
  gsl_vector_free( work );
  gsl_vector_free( b );
  gsl_vector_free( x );

  return EXIT_SUCCESS;
}

[-- Warning: decoded text below may be mangled, UTF-8 assumed --]
[-- Attachment #1.4: interp_test.c --]
[-- Type: text/x-csrc; name="interp_test.c", Size: 621 bytes --]

#include "find_interpolation_parameters.h"

int main( int argc, char * argv[] )
{
  const size_t dim = 2;
  const size_t num = 9;
  double point_data[] = { 
    1, 1,
    1, 2,
    1, 3,
    2, 1,
    2, 2,
    2, 3,
    3, 1,
    3, 2,
    3, 3 };
  
  double value_data[] = {
    2.00,
    2.11,
    3.21,
    3.97,
    3.76,
    4.01,
    3.40,
    5.47,
    4.35 };

  double parameters[3];

  find_interpolation_parameters_SVD( 
    dim, num, point_data, value_data, parameters );

  printf( "%f %f %f\n", parameters[0], parameters[1], parameters[2] );

  return EXIT_SUCCESS;
}

[-- Attachment #2: OpenPGP digital signature --]
[-- Type: application/pgp-signature, Size: 252 bytes --]