#include <stdio.h>
#include <math.h>
void product(int n, int m, double A[][100], double C[][100], double d[], double b[]);
void Cholesky(int n, double A[][100], double L[][100]);
void Chase(int m, double L[][100], double d[], double x[], double y[]);
int main(){
int n, m;
double A[100][100], C[100][100], d[100], b[100], L[100][100];
double x[100], y[100];
scanf("%d %d", &n, &m); //矩阵 A[][] n 行 m 列
for(int i = 0; i < n; i++){
for(int j = 0; j < m; j++){
scanf("%lf", &A[i][j]);
}
}
for(int i = 0; i < n; i++){
scanf("%lf", &b[i]);
}
product(n, m, A, C, d, b);
Cholesky(m, C, L);
Chase(m, L, d, x, y);
for(int i = 0; i < m; i++){
printf("%lf ", x[i]);
}
}
void product(int n, int m, double A[][100], double C[][100], double d[], double b[]){
//计算 C = A' * A
for (int i=0; i < n; i ++) {
for (int j=0; j < n; j ++) {
double ret = 0;
for(int k=0; k<n; k++) {
ret += A[k][i] * A[k][j];
}
C[i][j] += ret;
}
}
//计算 d = A' * b
for (int i=0; i<n; i ++) {
double ret = 0;
for(int j=0; j<n; j ++) {
ret += A[j][i] * b[j];
}
d[i] = ret;
}
}
void Cholesky(int n, double A[][100], double L[][100]){
// Cholesky 分解正定对称阵 A
for (int k =0; k<n; k ++) {
double temp = 0;
for (int j=0; j<k; j ++) {
temp += pow(L[k][j], 2);
}
L[k][k] = sqrt(A[k][k] - temp);
for (int i=k+1; i<n; i ++) {
double ret = 0;
for (int j=0; j<k; j ++) {
ret += L[i][j] * L[k][j];
}
L[i][k] = (A[i][k] - ret) / L[k][k];
}
}
// for (int i=0; i<n; i ++) {
// for (int j=0; j<n; j ++) {
// printf("%lf \t", L[i][j]);
// }
// printf("\n");
// }
}
void Chase(int n, double L[][100], double d[], double x[], double y[]){
//求解方程组 Ly = d
for (int i=0; i<n; i ++) {
double ret = 0;
for (int j=0; j<i; j ++) {
ret = L[i][j] * y[j];
}
y[i] = (d[i] - ret) / L[i][i];
}
// for (int i=0; i<n; i ++) {
// printf("%lf \t", y[i]);
// }
// printf("\n");
//求解方程组 L'x = y
for (int i=0; i<n; i ++) {
double ret = 0;
for (int j=0; j<i; j ++) {
ret = L[n-j-1][n-i-1] * x[n-j-1];
}
x[n-i-1] = (y[n-i-1] - ret) / L[n-i-1][n-i-1];
}
}
正则化方法
于 2023-04-12 12:08:46 首次发布