c++ - 使用 BLAS 和 LAPACKE 在 C++ 中使用 SVD 计算 Pseidoinverse
问题描述
我正在尝试实现矩阵的伪逆计算 A*,以便为具有 C++ 维度的方形 nxn 矩阵 A 求解 Ax=b。A* 的算术公式是通过 SVD 分解。
因此,首先我计算 SVD(A)=USV^T,然后计算 A*=VS U^T,其中 S是反对角线 S,其中其非零元素 si 在 S* 中变为 1/si。最后我计算解决方案 x=A*b
但是我没有得到正确的结果。我将 LAPACKE 接口用于 c++ 和 cblas 用于矩阵乘法。这是我的代码:
double a[n * n] = {2, -1, 2,1};
double b[n]={3,4};
double u[n * n], s[n],vt[n * n];
int lda = n, ldu = n, ldvt = n;
int info = LAPACKE_dgesdd(LAPACK_COL_MAJOR, 'A', n, n, a, lda, s,
u, ldu, vt, ldvt);
for (int i = 0; i < n; i++) {
s[i] = 1.0 / s[i];
}
const int a = 1;
const int c = 0;
double r1[n];
double r2[n];
double res[n];
//compute the first multiplication s*u^T
cblas_dgemm( CblasColMajor,CblasNoTrans, CblasTrans, n, n, n, a, u, ldvt, s, ldu, c, r1, n);
//compute the second multiplication v^T^T=vs*u^T
cblas_dgemm( CblasColMajor,CblasTrans, CblasNoTrans, n, n, n, a, vt, ldvt, r1, ldu, c, r2, n);
//now that we have the pseudoinverse A* solve by multiplying with b.
cblas_dgemm( CblasColMajor,CblasNoTrans, CblasNoTrans, n, 1, n, a, r2, ldvt, b, ldu, c, res, n);
在第二个 cblas_dgemm 之后,预计在 r2 中有 A* 伪逆。但是,在与 matlab pinv 比较后,我没有得到相同的结果。如果我打印 r2 结果给出:
0.25 0.50
0.25 0.50
但应该是
0.25 -0.50
0.25 0.50
解决方案
的参数表示S
LAPACKE_dgesdd()
SVD 分解中矩阵的奇异值。虽然它的长度为n
,但它没有描述向量,因为它表示对角矩阵。事实上, Su^T 的结果是一个大小为 的矩阵n*n
。
cblas_dscal()
可以在循环中应用该例程来计算涉及对角矩阵的矩阵乘积,尽管生成的 Su^t 仍然是转置的。看看在 fortran 中乘以对角矩阵的最佳方法是什么
以下代码可以通过g++ main.cpp -o main -llapacke -llapack -lgslcblas -lblas -lm -Wall
(或-lcblas`...)编译
#include <iostream>
#include <string>
#include <fstream>
#include <stdlib.h>
#include <stdio.h>
#include <math.h>
extern "C" {
#include <lapacke.h>
#include <cblas.h>
}
int main(int argc, char *argv[])
{
const int n=2;
double a[n * n] = {2, -1, 2,1};
double b[n]={3,4};
double u[n * n], s[n],vt[n * n];
int lda = n, ldu = n, ldvt = n;
//computing the SVD
int info = LAPACKE_dgesdd(LAPACK_COL_MAJOR, 'A', n, n, a, lda, s,
u, ldu, vt, ldvt);
if (info !=0){
std::cerr<<"Lapack error occured in dgesdd. error code :"<<info<<std::endl;
}
for (int i = 0; i < n; i++) {
s[i] = 1.0 / s[i];
}
const int aa = 1;
const int c = 0;
//double r1[n*n];
double r2[n*n];
double res[n];
//compute the first multiplication s*u^T
// here : s is not a vector : it is a diagonal matrix. The ouput must be of size n*n
//cblas_dgemm( CblasColMajor,CblasNoTrans, CblasTrans, n, n, n, aa, u, ldvt, s, ldu, c, r1, n);
for (int i = 0; i < n; i++) {
cblas_dscal(n,s[i],&u[i*n],1);
}
//compute the second multiplication v^T^T=vs*u^T
cblas_dgemm( CblasColMajor,CblasTrans, CblasTrans, n, n, n, aa, vt, ldvt, u, ldu, c, r2, n);
//now, r2 is the pseudoinverse of a.
//now that we have the pseudoinverse A* solve by multiplying with b.
cblas_dgemm( CblasColMajor,CblasNoTrans, CblasNoTrans, n, 1, n, aa, r2, ldvt, b, ldu, c, res, n);
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
std::cout<<r2[i*n+j]<<" ";
}
}
std::cout<<std::endl;
}
它打印预期的结果:
0.25 0.25 -0.5 0.5
推荐阅读
- c# - 如何降级 Autorest 扩展 - 特别是 C# 扩展?
- java - Android - 构建失败 - 任务':app:mergeReleaseResources'的执行失败
- r - 读取文件而不在 R 中提及其特定文件名
- android - 最后一个元素在 LazyColumn 中不可见。喷气背包组成
- database - 如何恢复oracle db脚本,使用分区时标准版可以使用
- javascript - 替换字符串中的数组值并根据值创建新字符串
- python - 在Python中将非常小的数字转换为更大的正值
- python - pandas value_counts() 返回不正确的计数
- visual-c++ - 无法使用 OutLook_Bar 样式更改 CMFCPropertySheet 中的选项卡位图
- reactjs - react-google-recaptcha 允许 CSP