LAPACK 的奇异值分解:大矩阵的问题
Singular Value Decomposition with LAPACK: problems with big matrices
我正在使用 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列的矩阵),执行编译脚本会报错:Segmentation fault (core dumped).
我试图在 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: 我使用的是 3.6.0 版本的 lapacke,以及 Ubuntu 16.04 LTS。
像这样在堆上分配矩阵解决了问题:
double* s = new double[n];
double* u = new double[ldu * m];
double* vt = new double[ldvt * n];
double* a = new double[lda * n];
我正在使用 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列的矩阵),执行编译脚本会报错:Segmentation fault (core dumped).
我试图在 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: 我使用的是 3.6.0 版本的 lapacke,以及 Ubuntu 16.04 LTS。
像这样在堆上分配矩阵解决了问题:
double* s = new double[n];
double* u = new double[ldu * m];
double* vt = new double[ldvt * n];
double* a = new double[lda * n];