A^tAx=A^tb
の解となる
AをQR分解すると
A^tAx=R^tQ^tQRx=R^tRx
A^tb=R^tQ^tb
となる
(R^t)^-1を両辺にかけて
Rx=Q^tbを求めれば良い
Rは三角なので簡単に求められる
(0, -2)
(1, 2,)
(2, 5)
(3, 11)
(4, 1)に対して3次式で近似してみる
最小二乗法のプログラムは↓
#include<iostream>
using namespace std;
/* m行n列 */
double** init_matrix(int m, int n) {
double** A = new double*[m];
for (int i = 0; i < m; i++) {
A[i] = new double[n]();
for (int j = 0; j < n; j++) {
cin >> A[i][j];
}
}
return A;
}
void print_matrix(double** A, int n, int m) {
cout << endl;
cout << endl;
for (int i = 0; i < m; i++) {
for (int j = 0; j < n; j++) {
cout << A[i][j] << " ";
}
cout << endl;
}
}
/*Aの左下三角のそれぞれの列ににハウスホルダー変換のuが入
int qr(double** A, int n, int m, double* y) {
double x, z, w;
for (int i = 0; i < n; i++) {
x = 0;
for (int j = i; j < m; j++) {
x += A[j][i] * A[j][i];
}
if (x < 0.00001) {
return 1;
}
x = sqrt(x);
if (A[i][i] > 0) {
y[i] = -x;
}
else {
y[i] = x;
}
z = 1 / sqrt(y[i] * (y[i] - A[i][i]));
/* ハウスホルダー変換のuをAの左下三角の列に格納 */
A[i][i] = z*(A[i][i] - y[i]);
for (int j = i + 1; j < m; j++) {
A[j][i] = A[j][i] * z;
}
/* Aの右下にハウスホルダー変換を行う計算
H=E-uu^tに対してHv=v-(u,v)uとなる．
まず，w=(u,v)計算する．
*/
for (int j = i + 1; j < n; j++) {
w = 0;
for (int k = i; k < m; k++) {
w += A[k][j] * A[k][i];
}
for (int k = i; k < m; k++) {
A[k][j] -= w*A[k][i];
}
}
}
return 0;
}
void lsm(double** A,double* y, double* b,int m,int n) {
/* Aの右下にハウスホルダー変換を行う計算
H=E-uu^tに対してHb=b-(u,b)uとなる．
まず，w=(u,v)計算する．
*/
double w;
for (int i = 0; i < n; i++) {
w = 0;
for (int j = i; j < m; j++) {
w += b[j] * A[j][i];
}
for (int j = i; j < m; j++) {
b[j] -= w*A[j][i];
}
}
cout << " b print" << endl;
for (int i = 0; i < n; i++) {
cout << b[i] << " ";
}
cout << endl;
/*下から代入して解く*/
for (int i = n - 1; i >= 0; i--) {
for (int j = i + 1; j < n; j++) {
b[i] -= b[j] * A[i][j];
}
b[i] = b[i] / y[i];
}
}
int main() {
int n, m;
cin >> m >> n;
double* y=new double[n]();
double** A = init_matrix(m, n);
double* b = new double[m];
for (int i = 0; i < m; i++) {
cin >> b[i];
}
print_matrix(A, n, m);
qr(A, n, m, y);
print_matrix(A, n, m);
cout << "y" << endl;
for (int i = 0; i < m; i++) {
cout << y[i] << " ";
}
cout << endl;
lsm(A, y, b, m, n);
cout << endl;
cout << "b" << endl;
for (int i = 0; i < n; i++) {
cout<< b[i]<<endl;
}
return 0;
}
結果
1 0 0 0
1 1 1 1
1 2 4 8
1 3 9 27
1 4 16 64
matend
1.203 -4.47214 -13.4164 -44.7214
0.371748 -1.05867 12.6491 48.6991
0.371748 0.184608 -1.33233 22.4499
0.371748 0.48331 -0.472748 -1.13192
0.371748 0.782012 0.0374653 0.847794
matend
y
-2.23607 3.16228 3.74166 3.79473 2.10535e-314
b print
-7.60263 4.74342 -6.68153 -4.74342
b
-1.67143
-2.10714
5.71429
-1.25
0 件のコメント:
コメントを投稿