用 GSL 求解超定方程组及矩阵的奇异值分解(SVD) 2

接上一篇。。。

下面我们将 SVD 相关的功能封装成一个类，以方便我们提取 $S$ 和 $V$ 的值。
另外，当我们一个 $A$ 有多组 $x$ 需要求解时，也只需要计算一次 SVD 分解，用下面的类能减少很多计算量。

头文件如下：

    #ifndef GSLSINGULARVALUEDECOMPOSITION_H
    #define GSLSINGULARVALUEDECOMPOSITION_H

    #include <gsl/gsl_matrix.h>
    #include <gsl/gsl_vector.h>
    #include <gsl/gsl_blas.h>
    #include <gsl/gsl_linalg.h>
    #include <gsl/gsl_errno.h>

    void linearSolve_SVD(const gsl_matrix * A, const gsl_vector * b, gsl_vector * x);

    class GslSVD
    {
    public:
        GslSVD();
        ~GslSVD();
        int SV_decomp(const gsl_matrix * A);
        int SV_decomp_mod(const gsl_matrix * A);
        int SV_decomp_jacobi (gsl_matrix * A);
        int SV_solve(const gsl_vector *b, gsl_vector *x);

        gsl_vector * getVectorS();
        gsl_matrix * getMatrixU();
        gsl_matrix * getMatrixV();

        int trimVectorS(double abseps);
    private:
        gsl_vector * S;
        gsl_matrix * U;
        gsl_matrix * V;

        void alloc_suv(int rows, int cols);
    };

    #endif // GSLSINGULARVALUEDECOMPOSITION_H

程序文件如下：

#include "gsl_SVD.h"

void linearSolve_SVD(const gsl_matrix * A, const gsl_vector * b, gsl_vector * x)
{
    int rows = A->size1;
    int cols = A->size2;
    gsl_vector * work = gsl_vector_alloc (cols);
    gsl_vector * S = gsl_vector_alloc (cols);
    gsl_matrix * U = gsl_matrix_alloc(rows, cols);;
    gsl_matrix * V = gsl_matrix_alloc(cols, cols);

    gsl_matrix_memcpy (U, A); // 为了不破坏 A 中原始的数据，这里全都拷贝到 U 中

    gsl_linalg_SV_decomp( U, V, S, work );
    gsl_linalg_SV_solve ( U, V, S, b, x );

    gsl_vector_free(work);
    gsl_vector_free(S);
    gsl_matrix_free(V);
    gsl_matrix_free(U);
}
int GslSVD::trimVectorS(double abseps)
{
    int count = 0;
    for(int i = 0; i < S->size; i++)
    {
        if(fabs(gsl_vector_get(S, i)) < abseps)
        {
            count ++;
            gsl_vector_set(S, i, 0);
        }
    }
    return count;
}

gsl_vector * GslSVD::getVectorS()
{
    if(S == NULL) return NULL;
    gsl_vector * s = gsl_vector_alloc(S->size);
    gsl_vector_memcpy(s, S);
    return s;
}

gsl_matrix * GslSVD::getMatrixU()
{
    if(U == NULL) return NULL;
    gsl_matrix * u = gsl_matrix_alloc(U->size1, U->size2);
    gsl_matrix_memcpy(u, U);
    return u;
}

gsl_matrix * GslSVD::getMatrixV()
{
    if(V == NULL) return NULL;
    gsl_matrix * v = gsl_matrix_alloc(V->size1, V->size2);
    gsl_matrix_memcpy(v, V);
    return v;
}

GslSVD::GslSVD()
{
    S = NULL;
    U = NULL;
    V = NULL;
}

void GslSVD::alloc_suv(int rows, int cols)
{
    if( S != NULL )
    {
        gsl_vector_free(S);
        gsl_matrix_free(U);
        gsl_matrix_free(V);
    }
    S = gsl_vector_alloc (cols);
    U = gsl_matrix_alloc(rows, cols);
    V = gsl_matrix_alloc(cols, cols);
}

int GslSVD::SV_decomp(const gsl_matrix * A)
{
    int rows = A->size1;
    int cols = A->size2;

    gsl_vector * work = gsl_vector_alloc (cols);

    alloc_suv(rows, cols);
    gsl_matrix_memcpy (U, A); // 为了不破坏 A 中原始的数据，这里全都拷贝到 U 中
    int ret = gsl_linalg_SV_decomp( U, V, S, work );

    gsl_vector_free(work);

    return ret;
}

int GslSVD::SV_decomp_mod(const gsl_matrix * A)
{
    int rows = A->size1;
    int cols = A->size2;

    gsl_vector * work = gsl_vector_alloc (cols);
    gsl_matrix *X = gsl_matrix_alloc(cols, cols);

    alloc_suv(rows, cols);
    gsl_matrix_memcpy (U, A); // 为了不破坏 A 中原始的数据，这里全都拷贝到 U 中
    int ret = gsl_linalg_SV_decomp_mod( U, X, V, S, work );

    gsl_matrix_free(X);
    gsl_vector_free(work);

    return ret;
}

int GslSVD::SV_decomp_jacobi (gsl_matrix * A)
{
    int rows = A->size1;
    int cols = A->size2;
    alloc_suv(rows, cols);
    gsl_matrix_memcpy (U, A); // 为了不破坏 A 中原始的数据，这里全都拷贝到 U 中
    int ret = gsl_linalg_SV_decomp_jacobi( U, V, S );
    return ret;
}

int GslSVD::SV_solve(const gsl_vector *b, gsl_vector *x)
{
    if(U != NULL)
    {
        return gsl_linalg_SV_solve (U, V, S, b, x);
    }
    return -1;
}

GslSVD::~GslSVD()
{
    if(S != NULL)
    {
        gsl_vector_free(S);
        gsl_matrix_free(V);
        gsl_matrix_free(U);
    }
}

下面是个简单的测试代码：

void test5()
{
    double a_data[] = {1, 2,
                       2, 4};
    gsl_matrix_view A = gsl_matrix_view_array (a_data, 2, 2);
    GslSVD svd;
    svd.SV_decomp(&A.matrix);

    puts("S = ");
    gsl_vector_fprintf (stdout, svd.getVectorS(), "%f");

    puts("\nV = ");
    gsl_matrix_fprintf (stdout, svd.getMatrixV(), "%f");

    double b_data[] = {3, 6};
    gsl_vector_view b = gsl_vector_view_array (b_data, 2);
    gsl_vector * x = gsl_vector_alloc (2);
    svd.SV_solve(&b.vector, x);

    puts("\nx = ");
    gsl_vector_fprintf (stdout, x, "%f");
}