A demonstration for calculating pseudo-inverse matrix using the GNU Scientific Library (GSL).

gsl_pinv.c

/* *****************************************************************************
A demonstration for calculating pseudo-inverse matrix using the GNU Scientific
Library (GSL).

Copyright (C) 2025 by Akira TAMAMORI

This program is free software; you can redistribute it and/or modify it under
the terms of the GNU General Public License as published by the Free Software
Foundation; either version 2 of the License, or (at your option) any later
version.

This program is distributed in the hope that it will be useful, but WITHOUT ANY
WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR A
PARTICULAR PURPOSE.  See the GNU General Public License for more details.

You should have received a copy of the GNU General Public License along with
this program; if not, write to the Free Software Foundation, Inc., 51 Franklin
Street, Fifth Floor, Boston, MA 02110-1301, USA.
***************************************************************************** */

/* *****************************************************************************
This program utilizes the GNU Scientific Library (GSL).  GSL is licensed under
the GNU General Public License version 3 or later.  See the GSL documentation or
source code for its copyright notice.
***************************************************************************** */

#include <gsl/gsl_linalg.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_vector.h>
#include <stdio.h>

double frobenius_norm(const gsl_matrix *A, const gsl_matrix *B) {
    /* Calculate Frobenius norm between two matrices.

    Args:
        A (const gsl_matrix *): The first matrix.
        B (const gsl_matrix *): The second matrix.
    */
    double sum = 0.0;
    for (size_t i = 0; i < A->size1; i++) {
        for (size_t j = 0; j < A->size2; j++) {
            double diff = gsl_matrix_get(A, i, j) - gsl_matrix_get(B, i, j);
            sum += diff * diff;
        }
    }
    return sqrt(sum);
}

void pseudo_inverse(const gsl_matrix *M, double tolerance_level,
                    gsl_matrix *M_pinv) {
    /* Calculate the pseudo-inverse for a given matrix.

    Args:
        M (const gsl_matrix *): The original matrix.
        tolerance_level (double): The tolerance level of singular value.
        M_pinv (gsl_matrix *): The pseudo-inverse of M.
    */
    int m = M->size1;
    int n = M->size2;
    gsl_matrix *U = gsl_matrix_alloc(m, n);
    gsl_matrix *V = gsl_matrix_alloc(n, n);
    gsl_vector *S = gsl_vector_alloc(n);
    gsl_vector *work = gsl_vector_alloc(n);
    gsl_vector *col_vector = gsl_vector_alloc(m);
    gsl_matrix_memcpy(U, M);
    gsl_linalg_SV_decomp(U, V, S, work);
    for (int i = 0; i < n; i++) {
        double s = gsl_vector_get(S, i);
        gsl_matrix_get_col(col_vector, U, i);
        if (s > tolerance_level) {
            gsl_vector_scale(col_vector, 1.0 / s);
        } else {
            gsl_vector_set_zero(col_vector);
        }
        gsl_matrix_set_col(U, i, col_vector);
    }
    gsl_blas_dgemm(CblasNoTrans, CblasTrans, 1.0, V, U, 0.0, M_pinv);
    gsl_matrix_free(U);
    gsl_matrix_free(V);
    gsl_vector_free(S);
    gsl_vector_free(work);
    gsl_vector_free(col_vector);
}

int main(void) {
    int m = 3;  // number of rows
    int n = 2;  // number of columns

    // Initialize matrix A
    gsl_matrix *A = gsl_matrix_alloc(m, n);
    gsl_matrix_set(A, 0, 0, 1.0);
    gsl_matrix_set(A, 0, 1, 2.0);
    gsl_matrix_set(A, 1, 0, 3.0);
    gsl_matrix_set(A, 1, 1, 4.0);
    gsl_matrix_set(A, 2, 0, 5.0);
    gsl_matrix_set(A, 2, 1, 6.0);

    // Calculate the pseudo-inverse A+
    gsl_matrix *A_pinv = gsl_matrix_alloc(n, m);
    double tol = 1e-10;  // Set an appropriate threshold
    pseudo_inverse(A, tol, A_pinv);

    // Calculate the reconstructed matrix: A * A+ * A = A
    gsl_matrix *A_reconstructed = gsl_matrix_alloc(m, n);
    gsl_matrix *temp_matrix = gsl_matrix_alloc(m, m);
    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, A, A_pinv, 0.0,
                   temp_matrix);  // A * A+
    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, temp_matrix, A, 0.0,
                   A_reconstructed);  // (A * A+) * A
    gsl_matrix_free(temp_matrix);

    // Calculate the reconstruction error: ||A - A_reconstructed||_F
    double reconstruction_error = frobenius_norm(A, A_reconstructed);

    // Display the results
    printf("Original Matrix A:\n");
    for (int i = 0; i < m; i++) {
        for (int j = 0; j < n; j++) {
            printf("%f ", gsl_matrix_get(A, i, j));
        }
        printf("\n");
    }

    printf("\nPseudo Inverse A+:\n");
    for (int i = 0; i < n; i++) {
        for (int j = 0; j < m; j++) {
            printf("%f ", gsl_matrix_get(A_pinv, i, j));
        }
        printf("\n");
    }

    // Display the reconstruction error
    printf("\nReconstruction Error: %e\n", reconstruction_error);

    // Free memory
    gsl_matrix_free(A);
    gsl_matrix_free(A_pinv);

    return 0;
}