* LinAlg SV Doc?
@ 2010-02-16 19:28 Kevin H. Hobbs
0 siblings, 0 replies; only message in thread
From: Kevin H. Hobbs @ 2010-02-16 19:28 UTC (permalink / raw)
To: gsl-discuss
[-- 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 --]
^ permalink raw reply [flat|nested] only message in thread
only message in thread, other threads:[~2010-02-16 19:28 UTC | newest]
Thread overview: (only message) (download: mbox.gz / follow: Atom feed)
-- links below jump to the message on this page --
2010-02-16 19:28 LinAlg SV Doc? Kevin H. Hobbs
This is a public inbox, see mirroring instructions
for how to clone and mirror all data and code used for this inbox;
as well as URLs for read-only IMAP folder(s) and NNTP newsgroup(s).