LLE.cpp

#include <stdio.h>
#include <gsl/gsl_matrix_double.h>
#include <gsl/gsl_vector_double.h>
#include <gsl/gsl_blas.h>
#include <gsl/gsl_linalg.h>
#include <gsl/gsl_cblas.h>
#include <gsl/gsl_eigen.h>
#include <gsl/gsl_math.h>
#include <gsl/gsl_sort_vector.h>
#include <gsl/gsl_sf.h>

/**
 * @param matrix The matrix waiting for transpose
 * @return Return transpose of input matrix
 */
gsl_matrix* viewDice (gsl_matrix* matrix) {
    int row = matrix->size1;
    int column = matrix->size2;
    gsl_matrix* output = gsl_matrix_alloc(column, row);
    for (int i=0; i<row; ++i){
        for (int j=0; j<column; ++j){
            double element = gsl_matrix_get(matrix, j, i);
            gsl_matrix_set(output, i, j, element);
        }
    }
    return output;
}

/**
 * @param matrix The matrix waiting for calculating invert matrix
 * @return Return the invert matrix of input matrix
 */
gsl_matrix* invertMatrix (gsl_matrix* matrix) {
    int row = matrix->size1;
    int column = matrix->size2;
    int s;
    if (row != column){
        return null;
    }
    gsl_matrix * inverse = gsl_matrix_alloc (row, column);
    gsl_permutation * perm = gsl_permutation_alloc (row);
    gsl_linalg_LU_decomp (matrix, perm, &s);
    gsl_linalg_LU_invert (matrix, perm, inverse);
    return inverse
}


void runLLE(gsl_matrix matrix, int P){

    /* P 代表 pivot 的个数 */
    const int numOfPivot = P;

    /* 最近邻点的个数，默认为 18 */
    const int nearestNeighbors = 18;

    /* 正规化参数 */
    const double regularizationParameter = 0.001;

    int numOfRow = m.size1;
    int numOfColumn = m.size2;

    /* 为权值矩阵申请内存空间 */
    gsl_matrix* weightMatrix = gsl_matrix_alloc(numOfRow, numOfColumn);
    for (int i=0; i<numOfRow; ++i){
        for (int j=; j<numOfColumn; ++j){
            gsl_matrix_set(weightMatrix, i, j, 0);
        }
    }

    /* 单位矩阵 */ 
    gsl_matrix* E = gsl_matrix_alloc(k, k);
    for (int i=0; i<k; ++i){
        gsl_matrix_set(E, i, i, 1);
    }

    for (int row=0; row<numOfRow; ++row){

        /* 步骤1， 找出 k 个最近的邻点 */
        gsl_matrix* diffMatrix = gsl_matrix_alloc(numOfRow, numOfColumn);

        /* 计算 [row, i] 和 [j, i] 之间的差，将差的平方存入 diffSquare 矩阵中 */
        for (int i=0; i<numOfColumn; ++i){
            double element = gsl_matrix_get(matrix, row, i);
            for (int j=0; j<numOfRow; ++j){
                double diff = gsl_matrix_get(matrix, j, i) - element;
                double diffSquare = diff * diff;
                gsl_matrix_set(diffMatrix, j, i, diffSquare);
            }
        }

        /* 计算 diffMatrix 每一列的元素之和 */
        gsl_vector* sumOfSquareVector = gsl_vector_alloc(numOfRow);
        gsl_vector* sortedIndex = gsl_vector_alloc(numOfRow);
        for (int i=0;i<numOfRow;++i){
            gsl_vector_set(sortedIndex, i, i);
        }

        for (int i=0; i<numOfRow; ++i){
            double element = 0;
            for (int j=0; j<numOfColumn; ++j){
                element = element + gsl_matrix_get(diffMatrix, i, j);
            }
            gsl_vector_set(sumOfSquareVector, i, element);
        }

        gsl_sort_vector2(sumOfSquareVector, sortedIndex);

        gsl_matrix* sortedMatrix = gsl_matrix_alloc(nearestNeighbors, numOfColumn);
        for (int i=1; i<=nearestNeighbors; ++i){
            int rowTemp = gsl_vector_get(sortedIndex, i);
            gsl_matrix_get_row(sortedMatrix, matrix, rowTemp);
        }

        /* 步骤2: 重建矩阵 */
        gsl_matrix* sortedMatrixDice = viewDice(sortedMatrix);
        gsl_matrix* Q = gsl_matrix_alloc(nearestNeighbors, nearestNeighbors);
        
        /* 下面的语句是计算 Q = sortedMatrix x sortedMatrixDice */
        gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, sortedMatrix, sortedMatrixDice, 0.0, Q);

        double matrixTrace = 0;
        for (int i=0; i<nearestNeighbors; ++i){
            matrixTrace = matrixTrace + gsl_matrix_get(Q, i, i);
        }

        double base = regulariztionParameter * matrixTrace;
        for (int Qi=0; Qi<nearestNeighbors; ++Qi){
            double element = gsl_matrix_get(Q, Qi, Qi);
            gsl_matrix_set(Q, Qi, Qi, element + base);
        }
        gsl_matrix* invertQ = invertMatrix(Q);
        gsl_matrix* wTemp = gsl_matrix_alloc(k, k);
        gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, invertQ, E, 0.0, wTemp);

        gsl_vector* w = gsl_vector_alloc(k);
        gsl_matrix_get_col(w, wTemp, 0);
        double wSum = 0;
        for (int i=0;i<k;++i){
            wSum += gsl_vector_get(w, i);
        }

        for (int i=0;i<w.size;++i){
            double temp = gsl_vector_get(w, i) / wSum;
            gsl_vector_set(w, i, temp);
        }

        for (int i=0;i<k;++i){
            double temp = gsl_vector_get(w, i);
            int index = (int) gsl_vector_get(sortedIndex, i+1);
            gsl_matrix_set(weightMatrix, index, row, temp);
        }

        gsl_matrix_free(diffMatrix);
        gsl_matrix_free(sortedMatrix);
        gsl_matrix_free(sortedMatrixDice);
        gsl_matrix_free(Q);
        gsl_vector_free(sumOfSquareVector);
        gsl_vector_free(sortedIndex);
    }

    /* 步骤3：映射到低维空间中 */
    for (int i=0;i<weightMatrix.size1;++i){
        double newVal = gsl_matrix_get)weightMatrix, i, i) - 1;
        gsl_matrix_set(weightMatrix, i, i, newVal);
    }

    gsl_matrix* weightMatrixDice = viewDice(weightMatrix);
    gsl_matrix* M = gsl_matrix_alloc(weightMatrix.size1, weightMatrix.size2);
    gsl_blas_dgemm(CblasNoTrans, CblasNoTrans, 1.0, weightMatrix, weightMatrixDice, 0.0, M);

    for (int i=0;i<M.size2;++i){
        double element = gsl_matrix_get(M, i, i) + 0.1;
        gsl_matrix_set(M, i, i, element);
    }
    
    gsl_matrix * output = gsl_matrix_alloc (numOfRow, numOfPivot);
    
    /* eval 存储 N 个特征值 */
    /* evec 为由 M 矩阵的 N 的特征值向量构成的 N 阶方阵 */
    gsl_vector * eval = gsl_vector_alloc(numOfRow);
    gsl_matrix * evec = gsl_matrix_alloc(numOfRow, numOfRow);
    /* 申请工作空间 */
    gsl_eigen_symmv_workspace *w = gsl_eigen_symmv_alloc(numOfRow); 
    /* eval，evec 赋值 */
    gsl_eigen_symmv (matrix, eval, evec, w);

    /* eval 升序排序，对应 eval 的顺序给 evec 重新赋值 */
    gsl_eigen_symmv_sort(eval, evec, GSL_EIGEN_SORT_ABS_ASC);

    /* 需要转置矩阵evec使得行列正确对应 */
    gsl_matrix * evec2 = viewDice(evec);
    int i, j;
    for(i=0;i<numOfRow;i++){
        for(j=1;j<numOfPivot+1;j++){
            double temp = gsl_matrix_get(evec2, i, j);
            gsl_matrix_set(output, i, j-1, temp);
        }
    }

    gsl_eigen_symmv_free(w);
    gsl_vector_free(eval);
    gsl_matrix_free(evec);
    gsl_matrix_free(weightMatrix);
    gsl_matrix_free(E);

    return output;
}