Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Compute the (Moore-Penrose) pseudo-inverse of a libgsl matrix in plain C.
/**
* Compute the (Moore-Penrose) pseudo-inverse of a libgsl matrix in plain C.
*
* Compile uding:
*
* gcc moore_penrose_pseudoinverse.c -lgsl -lblas
*
* Dependencies:
* - libgsl (GNU Scientific Library)
* - libblas (Basic Linear Algebra Subprograms)
*
* Charl Linssen <charl@itfromb.it>
* Feb 2016
* PUBLIC DOMAIN
**/
#include <stdio.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
typedef double realtype;
#define max(a,b) ((a) > (b) ? (a) : (b))
#define min(a,b) ((a) < (b) ? (a) : (b))
typedef unsigned int bool;
#define false 0
#define true 1
void print_matrix(const gsl_matrix *m) {
size_t i, j;
for (i = 0; i < m->size1; i++) {
for (j = 0; j < m->size2; j++) {
printf("%f\t", gsl_matrix_get(m, i, j));
}
printf("\n");
}
}
void print_vector (const gsl_vector * v) {
size_t i;
for (i = 0; i < v->size; i++) {
printf("%f\t", gsl_vector_get (v, i));
}
}
/**
* Compute the (Moore-Penrose) pseudo-inverse of a matrix.
*
* If the singular value decomposition (SVD) of A = UΣVᵀ then the pseudoinverse A⁻¹ = VΣ⁻¹Uᵀ, where ᵀ indicates transpose and Σ⁻¹ is obtained by taking the reciprocal of each nonzero element on the diagonal, leaving zeros in place. Elements on the diagonal smaller than ``rcond`` times the largest singular value are considered zero.
*
* @parameter A Input matrix. **WARNING**: the input matrix ``A`` is destroyed. However, it is still the responsibility of the caller to free it.
* @parameter rcond A real number specifying the singular value threshold for inclusion. NumPy default for ``rcond`` is 1E-15.
*
* @returns A_pinv Matrix containing the result. ``A_pinv`` is allocated in this function and it is the responsibility of the caller to free it.
**/
gsl_matrix* moore_penrose_pinv(gsl_matrix *A, const realtype rcond) {
gsl_matrix *V, *Sigma_pinv, *U, *A_pinv;
gsl_matrix *_tmp_mat = NULL;
gsl_vector *_tmp_vec;
gsl_vector *u;
realtype x, cutoff;
size_t i, j;
unsigned int n = A->size1;
unsigned int m = A->size2;
bool was_swapped = false;
if (m > n) {
/* libgsl SVD can only handle the case m <= n - transpose matrix */
was_swapped = true;
_tmp_mat = gsl_matrix_alloc(m, n);
gsl_matrix_transpose_memcpy(_tmp_mat, A);
A = _tmp_mat;
i = m;
m = n;
n = i;
}
/* do SVD */
V = gsl_matrix_alloc(m, m);
u = gsl_vector_alloc(m);
_tmp_vec = gsl_vector_alloc(m);
gsl_linalg_SV_decomp(A, V, u, _tmp_vec);
gsl_vector_free(_tmp_vec);
/* compute Σ⁻¹ */
Sigma_pinv = gsl_matrix_alloc(m, n);
gsl_matrix_set_zero(Sigma_pinv);
cutoff = rcond * gsl_vector_max(u);
for (i = 0; i < m; ++i) {
if (gsl_vector_get(u, i) > cutoff) {
x = 1. / gsl_vector_get(u, i);
}
else {
x = 0.;
}
gsl_matrix_set(Sigma_pinv, i, i, x);
}
/* libgsl SVD yields "thin" SVD - pad to full matrix by adding zeros */
U = gsl_matrix_alloc(n, n);
gsl_matrix_set_zero(U);
for (i = 0; i < n; ++i) {
for (j = 0; j < m; ++j) {
gsl_matrix_set(U, i, j, gsl_matrix_get(A, i, j));
}
}
if (_tmp_mat != NULL) {
gsl_matrix_free(_tmp_mat);
}
/* two dot products to obtain pseudoinverse */
_tmp_mat = gsl_matrix_alloc(m, n);
gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1., V, Sigma_pinv, 0., _tmp_mat);
if (was_swapped) {
A_pinv = gsl_matrix_alloc(n, m);
gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1., U, _tmp_mat, 0., A_pinv);
}
else {
A_pinv = gsl_matrix_alloc(m, n);
gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1., _tmp_mat, U, 0., A_pinv);
}
gsl_matrix_free(_tmp_mat);
gsl_matrix_free(U);
gsl_matrix_free(Sigma_pinv);
gsl_vector_free(u);
gsl_matrix_free(V);
return A_pinv;
}
int main() {
const unsigned int N = 2;
const unsigned int M = 3;
const realtype rcond = 1E-15;
gsl_matrix *A = gsl_matrix_alloc(N, M);
gsl_matrix *A_pinv;
gsl_matrix_set(A, 0, 0, 1.);
gsl_matrix_set(A, 0, 1, 3.);
gsl_matrix_set(A, 0, 2, 5.);
gsl_matrix_set(A, 1, 0, 2.);
gsl_matrix_set(A, 1, 1, 4.);
gsl_matrix_set(A, 1, 2, 6.);
printf("A matrix:\n");
print_matrix(A);
A_pinv = moore_penrose_pinv(A, rcond);
printf("\nPseudoinverse of A:\n");
print_matrix(A_pinv);
gsl_matrix_free(A);
gsl_matrix_free(A_pinv);
return 0;
}
@jehfuger

This comment has been minimized.

Copy link

@jehfuger jehfuger commented Aug 27, 2017

I am trying to use the code above in cygwin and despite the fact that I have no error associated, I can`t get any answer.
Can someone help me?

@GhostNaN

This comment has been minimized.

Copy link

@GhostNaN GhostNaN commented Sep 28, 2020

You forgot:
#include <gsl/gsl_linalg.h>
For:
gsl_linalg_SV_decomp(A, V, u, _tmp_vec);
https://gist.github.com/turingbirds/5e99656e08dbe1324c99#file-moore_penrose_pseudoinverse-c-L92

@turingbirds

This comment has been minimized.

Copy link
Owner Author

@turingbirds turingbirds commented Sep 28, 2020

Thank you @GhostNaN! Gist has been updated.

@GhostNaN

This comment has been minimized.

Copy link

@GhostNaN GhostNaN commented Sep 29, 2020

No problem @turingbirds! It works great for getting my Savitzky Golay coefficients.

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment