Вычисление Pseidoinverse с SVD в C ++ с использованием BLAS и LAPACKE - PullRequest
Новая
       104

Вычисление Pseidoinverse с SVD в C ++ с использованием BLAS и LAPACKE

1 голос
/ 09 апреля 2019

Я пытаюсь реализовать псевдообратное вычисление A * матрицы для решения Ax = b для квадратной nxn матрицы A с размерами в C ++. Арифметическая формула для A * через разложение SVD.

Итак, сначала я вычисляю SVD (A) = USV ^ T, а затем A * = VS U ^ T, где S - обратная диагональ S, где ее ненулевой элемент si становится 1 / si в S *. Наконец я вычисляю решение 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 ожидается, что A * будет псевдообратным в r2. Однако после сравнения с Matlab Pinv я не получаю тот же результат. Если я печатаю r2, результат дает: 

 0.25   0.50
 0.25   0.50

но это должно быть 

0.25   -0.50
0.25   0.50

1 Ответ

0 голосов
/ 20 апреля 2019

Аргумент S из LAPACKE_dgesdd() представляет особые значения матрицы в разложении SVD . Хотя он имеет длину n, он не изображает вектор, поскольку представляет диагональную матрицу. Действительно, результатом S.u ^ T является матрица размера n*n.

Подпрограмма cblas_dscal() может быть применена в цикле для вычисления матричного произведения, включающего диагональную матрицу, хотя результирующий S.u ^ t все еще транспонирован. См. , как лучше всего умножить диагональную матрицу в Фортране

Следующий код может быть скомпилирован с помощью 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
...