'Moore-Penrose Pseudoinverse of matrix using Singular Value Decomposition
I am trying to compute the pseudoinverse of a matrix in C code. Specifically, I am trying to compute the Moore-Penrose pseudoinverse using singular value decomposition (SVD). I know there are some libraries out there for SVD but I am trying not to use any external libraries such as LAPACK or OpenCV. I found the following source code online and have been trying to follow it:
http://www.mymathlib.com/c_source/matrices/linearsystems/singular_value_decomposition.c
Sorry, it seems pretty long but a lot of it is good documentation. I am particularly interested in "After performing the singular value decomposition for A, call Singular_Value_Decomposition_Inverse to calculate the pseudo-inverse of A." Where I am struggling (I think mostly due to my inexperience with c) is how to go about running this source code.
My current main function:
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", )
}
I am not sure if I am on the right track so any guidance would be great. I am trying to understand what I am doing so please refrain from posting a completed solution. Just looking for some guidance :)
I do not understand what is meant by the following:
// 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. //
"After performing the singular value decomposition for A" Don't I have to call the Singular_Value_Decomposition first or are they referring to Singular_Value_Decomposition_Solve and it is a typo?
The problem I am currently facing is I'm not sure how to use the source code provided to generate an output. I understand the logic to calculate pinv is completed but I am not sure how to actually put it all together. I have started off with creating a main function to initialize some variables and my input matrix. After my initialization I am not sure what to do next
Solution 1:[1]
#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 available\n"); 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 converge\n");
else { printf(" The singular value decomposition of A is \n");
//printf("%d", err);
/*
for(int i=0;i<4;i++){
printf("%f\n", 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;
}
output:
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 . . .
Sources
This article follows the attribution requirements of Stack Overflow and is licensed under CC BY-SA 3.0.
Source: Stack Overflow
| Solution | Source |
|---|---|
| Solution 1 | kyazgan |