/*--------------------------------------------------------------------------------
正規方程式
normal equation
normal_eq.c
--------------------------------------------------------------------------------*/
#include "least_sq.h"
int least_sq_with_normal_eq( double y[ROW_NUM], double X[ROW_NUM][COL_NUM], double A[COL_NUM] )
{
int i;
int j;
int k;
double temp;
double XT[COL_NUM][ROW_NUM]; // Xの転置行列, transposed matrix of X
double temp_X[COL_NUM][COL_NUM];
double temp_XT[COL_NUM][ROW_NUM];
// 転置行列の作成, creating a transposed matrix
for( i = 0; i < ROW_NUM; i++ ){
for( j = 0; j < COL_NUM; j++ ){
XT[j][i] = X[i][j];
}
}
// 正規方程式による最小二乗法, Least squares method with normal equations
// Y = X * A
// A = (XT * X)^-1 * XT * Y
// temp_X = XT * X
for( i = 0; i < COL_NUM; i++ ){
for( j = 0; j < COL_NUM; j++ ){
temp = 0.0;
for( k = 0; k < ROW_NUM; k++ ){
temp += XT[i][k] * X[k][j];
}
temp_X[i][j] = temp;
}
}
// XT * Xの逆行列の計算, calculation of inverse matrix of XT * X
// 元々、temp_XにはXT * Xが入っていて、inv_matrixで(XT * X)^-1が返ってくる。
// Originally, temp_X contains XT * X and inv_matrix returns (XT * X)^-1.
if ( inv_matrix( temp_X ) == -1 ){
printf("# inverse matrix cannot be obtained.");
return -1;
}
// (XT * X)^-1 * XT
// temp_X = (XT * X)^-1
// temp_X * XT
for( i = 0; i < COL_NUM; i++ ){
for( j = 0; j < ROW_NUM; j++ ){
temp = 0.0;
for( k = 0; k < COL_NUM; k++ ){
temp += temp_X[i][k] * XT[k][j];
}
temp_XT[i][j] = temp;
}
}
// (XT * X)^-1 * XT * Y
// temp_XT * Y
for( i = 0; i < COL_NUM; i++ ){
temp = 0.0;
for( j = 0; j < ROW_NUM; j++ ){
temp += temp_XT[i][j] * y[j];
}
A[i] = temp;
}
return 0;
}