LAPACK 的奇异值分解：大矩阵的问题

Question

我正在使用 LAPACK 的 C 接口来计算矩阵的奇异值分解 (SVD)。为此，我正在使用例程

dgesvd_

。

我创建了一个简单的 C++ 脚本，它创建一个随机矩阵（具有

M

行和

N

列），并计算其 SVD。该脚本的代码如下：

#include <iostream>
#include <random>
#include "lapacke.h"

#define M 150
#define N 3
#define MEAN 0
#define STD 0.25


int main(int argc, char *argv[])
{

    std::default_random_engine generator;
    std::normal_distribution<double> distribution(MEAN, STD);

    int m = M, n = N, lda = M, ldu = M, ldvt = N, info, lwork;
    double wkopt;
    double* work;
    double s[n], u[ldu * m], vt[ldvt * n], a[lda * n];

    for(unsigned int k = 0; k < (lda * n); k++){
        a[k] = distribution(generator);
    }

    lwork = -1;
    dgesvd_("All", "All", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, &wkopt, &lwork, &info);
    lwork = (int)wkopt;
    work = (double*) malloc(lwork * sizeof(double));
    dgesvd_("All", "All", &m, &n, a, &lda, s, u, &ldu, vt, &ldvt, work, &lwork, &info);  

    if(info > 0){
        std::cerr<< "The algorithm computing SVD failed to converge\n"     ;
        exit(1);
    }

    free(work);

    std::cout<<"Singular values:\n";
    for(unsigned int k = 0; k < n; k++){
        std::cout<<s[k]<<" ";
    }
    std::cout<<"\n";
    return(0);
}

150 行 3 列的矩阵似乎可以正确计算 SVD。但是，当矩阵的行数较多时（例如，矩阵有 1500 行 3 列），编译脚本的执行会引发此错误：分段错误（核心已转储）。

我尝试在 Python 脚本中做类似的事情：

import numpy as np
M = 1500
N = 3
MEAN = np.zeros(3)
STD = 0.25
result = np.linalg.svd(points)

当我使用 Python 时，奇异值似乎可以正确计算，没有任何错误，尽管我使用的 NumPy 方法是使用 LAPACK 例程执行的

dgesdd_

（我尝试了这个，但出现了相同的错误）。

有谁知道为什么会发生这个错误？任何解决该问题的帮助将不胜感激。谢谢。

Pd：我使用的是 lapacke 3.6.0 版本，以及 Ubuntu 16.04 LTS。

Answer 1

像这样在堆上分配矩阵可以解决问题：

double* s = new double[n];
double* u = new double[ldu * m];
double* vt = new double[ldvt * n];
double* a = new double[lda * n];

Answer 2

打开

lapacke.h

，您可以看到负责工作计算的包装器。这很方便，因为这些可以在版本之间发生变化。开销可以忽略不计。

lapack_int LAPACKE_sgesdd( int matrix_layout, char jobz, lapack_int m,
                           lapack_int n, float* a, lapack_int lda, float* s,
                           float* u, lapack_int ldu, float* vt,
                           lapack_int ldvt );

第一个参数是 LAPACK_COL_MAJOR 或 LAPACK_ROW_MAJOR，具体取决于您的语言/平台。 此处显示了一个示例。