正则化方法

#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];
	}
}