我试图在C代码中计算矩阵的伪逆。具体来说，我试图使用奇异值分解(SVD)来计算Moore-Penrose伪逆。我知道有一些用于SVD的库，但我尽量不使用任何外部库，如LAPACK或OpenCV。我在网上找到了以下源代码，并一直试图遵循它:
http://www.mymathlib.com/c_source/matrices/linearsystems/singular_value_decomposition.c
对不起，它看起来很长，但其中很多都是好的文档。我特别感兴趣的是"在对A执行奇异值分解后，调用Singular_Value_Decomposition_Inverse来计算A的伪逆"。在我挣扎的地方(我认为主要是由于我对c缺乏经验)是如何去运行这个源代码。
我当前的主要功能:
int main(){
double A[4][4] = {{0, 0, 0, 0,},
{0, 2, 1, 2},
{2, 1, 0, 1},
{2, 0, 1, 4}};
#define M 4 //
#define N 4 //
double U[M][N]; //
double V[N][N]; //
double D[N]; //
double Astar[N][M]; //
double tolerance;
// here add whatever function need to be run
Singular_Value_Decomposition_Inverse(double* U, double* D, double* V,
double tolerance, int nrows, int ncols, double *Astar);
int i,j,k;
double *pu, *pv, *pa;
double dum;
dum = DBL_EPSILON * D[0] * (double) ncols;
if (tolerance < dum) tolerance = dum;
for ( i = 0, pv = V, pa = Astar; i < ncols; i++, pv += ncols)
for ( j = 0, pu = U; j < nrows; j++, pa++)
for (k = 0, *pa = 0.0; k < ncols; k++, pu++)
if (D[k] > tolerance) *pa += *(pv + k) * *pu / D[k];
printf(" The pseudo-inverse of A = UDV' is n", )
}
我不确定我是否在正确的轨道上，所以任何指导将是伟大的。我试图了解我在做什么，所以请避免张贴一个完整的解决方案。只是在寻找一些指导:)
我不明白下面这句话是什么意思:
// After performing the singular value decomposition for A, call //
// Singular_Value_Decomposition to solve the equation Ax = B or call //
// Singular_Value_Decomposition_Inverse to calculate the pseudo-inverse //
// of A. //
"对A"我不需要先调用Singular_Value_Decomposition吗?或者他们引用的是Singular_Value_Decomposition_Solve，这是一个错字?
我目前面临的问题是我不确定如何使用提供的源代码来生成输出。我理解计算pinv的逻辑是完成的，但我不确定如何实际地把它放在一起。我已经开始创建一个main函数来初始化一些变量和我的输入矩阵。初始化后，我不确定下一步该做什么
#include <stdio.h> // add
#include<stdlib.h> // add
#include <string.h> // required for memcpy()
#include <float.h> // required for DBL_EPSILON
#include <math.h> // required for fabs(), sqrt();
#define MAX_ITERATION_COUNT 30 // Maximum number of iterations
...
...
...
int main(){
int M=4;
int N=4;
double A[4][4] = {{0, 0, 0, 0},
{0, 2, 1, 2},
{2, 1, 0, 1},
{2, 0, 1, 4} };
double U[M][N];
double V[N][N];
double singular_values[N];
double* dummy_array;
// (your code to initialize the matrix A)
dummy_array = (double*) malloc(N * sizeof(double));
if (dummy_array == NULL) {printf(" No memory availablen"); exit(0); }
double err = Singular_Value_Decomposition((double*) A, M, N, (double*) U, singular_values, (double*) V, dummy_array); //
free(dummy_array);
if (err < 0) printf(" Failed to convergen");
else { printf(" The singular value decomposition of A is n");
//printf("%d", err);
/*
for(int i=0;i<4;i++){
printf("%fn", singular_values[i]);
}
*/
/*
for(int i=0;i<M;i++){
for(int j=0;j<M;j++){
//printf("t%0.3f" , A[i][j]);
printf("t%0.3f" , U[i][j]);
//printf("t%0.3f" , V[i][j]);
}
//printf("n");
//printf("n");
printf("n");
}
*/
}
double D[4]={5.408, 2.053, 1.59, 0};
double Astar[N][M];
double tolerance;
/*
double D[4][4] = {{5.408, 0, 0, 0},
{0, 2.053, 0, 0},
{0, 0, 0, 1.590},
{0, 0, 0, 0} };
*/
Singular_Value_Decomposition_Inverse((double*) U, D, (double*) V,
tolerance, M, N, (double*) Astar);
printf(" The pseudo-inverse of A = UDV' is n");
for(int i=0;i<M;i++){
for(int j=0;j<M;j++){
printf("t%0.6f" , Astar[i][j]);
}
printf("n");
}
return 0;
}
输出:
The pseudo-inverse of A = UDV' is
0.000000 -0.217970 0.461692 0.012779
0.000000 0.359000 0.269341 -0.256465
-0.000000 0.128219 -0.153902 0.051306
-0.000000 0.076937 -0.192391 0.230824
--------------------------------
Process exited after 0.06935 seconds with return value 0
Press any key to continue . . .