/*--------------------------------------------------------------------------------
QR分解を用いた最小二乗法
least square method using QR decomposition
main.c
--------------------------------------------------------------------------------*/
#include "main.h"
int main(void)
{
int i;
double A[COL_NUM]; // y = X * A
double y[ROW_NUM]; // y = X * A
double X[ROW_NUM][COL_NUM]; // y = X * A
// XをQR分解によってQとRに分解
// decompose X into Q and R by QR decomposition
double Q[ROW_NUM][ROW_NUM]; // X = Q * R
double R[ROW_NUM][COL_NUM]; // X = Q * R
/*----------------------------------------------------------------------*/
// y = X * A
// correct values of A
// A[0] = 1.0
// A[1] = 2.0
y[0] = 3.0;
y[1] = (double)DELTA;
y[2] = 2.0 * DELTA;
X[0][0] = 1.0; X[0][1] = 1.0;
X[1][0] = (double)DELTA; X[1][1] = 0.0;
X[2][0] = 0.0; X[2][1] = (double)DELTA;
/*----------------------------------------------------------------------*/
// QR分解の実行
// perform QR decomposition
householder_QR_decomp( X, Q, R );
// QR分解の結果の表示
// show X, Q and R (Q and R are the result of QR decomposition)
// show_XQR( X, Q, R );
// 最小二乗法（R*A=QT*yを用いてAの値を計算。QTはQの転置行列。）
// least square method (calculate the value of A using R*A=QT*y, where QT is the transpose matrix of Q)
solve_A_using_QR( y, Q, R, A );
printf("----------------------------------------------------\n");
printf(" number of rows: %d\n", (int)ROW_NUM);
printf(" number of columns: %d\n", (int)COL_NUM );
printf(" DELTA: %g\n", (double)DELTA );
printf("----------------------------------------------------\n");
for( i = 0; i < COL_NUM; i++ ){
printf(" A%d: %.10lf\n", i, A[i] );
}
printf("----------------------------------------------------\n");
return 0;
}
/*----------------------------------------------------------------------*/
void show_XQR( double X[ROW_NUM][COL_NUM], double Q[ROW_NUM][ROW_NUM], double R[ROW_NUM][COL_NUM] )
{
int i;
int j;
int k;
double s;
double t;
double temp;
// X[ROW_NUM][COL_NUM]
printf("Matrix X:\n");
for (i = 0; i < ROW_NUM; i++) {
for (j = 0; j < COL_NUM; j++) {
printf("%lf\t", X[i][j]);
// printf("%.16g\t", X[i][j]);
}
printf("\n");
}
printf("\n");
// Q[ROW_NUM][ROW_NUM]
printf("Matrix Q:\n");
for (i = 0; i < ROW_NUM; i++) {
for (j = 0; j < ROW_NUM; j++) {
printf("%lf\t", Q[i][j]);
// printf("%.16g\t", Q[i][j]);
}
printf("\n");
}
printf("\n");
// R[ROW_NUM][COL_NUM]
printf("Matrix R:\n");
for (i = 0; i < ROW_NUM; i++) {
for (j = 0; j < COL_NUM; j++) {
printf("%lf\t", R[i][j]);
// printf("%.16g\t", R[i][j]);
}
printf("\n");
}
printf("\n");
// Q * R
printf("Matrix Q * R:\n");
t = 0.0;
for (int i = 0; i < ROW_NUM; i++) {
for (int j = 0; j < COL_NUM; j++) {
s = 0.0;
for (k = 0; k < ROW_NUM; k++) {
// for (k = 0; k <= j; k++) {
s += Q[i][k] * R[k][j];
}
printf("%lf\t", s);
// printf("%.16g\t", s);
s -= X[i][j];
t += s * s;
}
printf("\n");
}
temp = t/(double)(ROW_NUM * COL_NUM);
// printf("# 二乗平均平方根誤差: %g\n", sqrt(temp) );
printf("# root-mean-square error: %g\n", sqrt(temp) );
printf("\n");
}