/usr/include/vmmlib/linear_least_squares.hpp is in libvmmlib-dev 1.0-2.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 | #ifndef __VMML__LINEAR_LEAST_SQUARES__HPP__
#define __VMML__LINEAR_LEAST_SQUARES__HPP__
#include <vmmlib/vector3.h>
#include "matrix_mxn.hpp"
namespace vmml
{
// beta_hat = ( Xt * X )^-1 * Xt * y;
// beta_hat -> best-fit 'solutions' to a,b = [ a,b ]transposed.
// TODO FIXME what about the constant???
struct linear_least_squares
{
template< size_t number_of_data_points, size_t number_of_unknowns, typename float_t >
vmml::matrix_mxn< number_of_unknowns, 1, float_t >
solve( std::vector< Vector3< float_t > >& data_points, float_t tolerance = 1e-9 )
{
assert( data_points.size() >= number_of_data_points );
matrix_mxn< number_of_data_points, number_of_unknowns, float_t > X;
matrix_mxn< number_of_data_points, 1, float_t > y;
for( size_t nb_index = 0; nb_index < number_of_data_points; ++nb_index )
{
float_t x = data_points[ nb_index ].x;
for( size_t index = 0; index < number_of_unknowns; ++index )
{
X[ nb_index ][ index ] = x;
x *= x;
}
y[ nb_index ][ 0 ] = data_points[ nb_index ].y;
}
std::cout << X << std::endl;
vmml::matrix_mxn< number_of_unknowns, number_of_data_points > Xt;
X.transposeTo( Xt );
std::cout << Xt << std::endl;
vmml::matrix_mxm< number_of_unknowns > XtX;
XtX.mul( Xt, X );
vmml::matrix_mxm< number_of_unknowns > XtX_inverse;
XtX.computeInverse( XtX_inverse, tolerance );
vmml::matrix_mxn< number_of_unknowns, number_of_data_points > XtX_inv_mul_Xt;
XtX_inv_mul_Xt.mul( XtX_inverse, Xt );
vmml::matrix_mxn< number_of_unknowns, 1 > beta_hat;
beta_hat.mul( XtX_inv_mul_Xt, y );
std::cout << beta_hat << std::endl;
return beta_hat;
}
template< size_t number_of_data_points, size_t number_of_unknowns, typename float_t >
vmml::matrix_mxn< number_of_unknowns, 1, float_t >
solve_using_qr_decomposition( std::vector< Vector3< float_t > >& data_points, float_t tolerance = 1e-9 )
{
assert( data_points.size() >= number_of_data_points );
matrix_mxn< number_of_data_points, number_of_unknowns, float_t > X;
matrix_mxn< number_of_data_points, 1, float_t > y;
// r = y - X * Beta
// X = Q*R
// Qt * r = Qt * y - ( Qt * Q ) R * Beta =
//
// | ( Qt * y )_n - R_n * Beta | | U |
// | ( Qt * y )_m-n | = | L |
//
// S = rt * Q * Qt * r = rt * r
//
// S = Ut * U + Lt * L
//
// => R_n * BetaHat = ( Qt * y )_n
//
vmml::matrix_mxn< number_of_unknowns, 1 > beta_hat;
return beta_hat;
}
template< size_t number_of_data_points, size_t number_of_unknowns, typename float_t >
vmml::matrix_mxn< number_of_unknowns, 1, float_t >
solve_using_svd( std::vector< Vector3< float_t > >& data_points, float_t tolerance = 1e-9 )
{
}
};
} // namespace vmml
#endif
|