Matrix Inverse using AFL & AQL


#1

Good morning,

I would like to use SciDB (14.12) to pseudo-invert matrices. I came up with the instructions below for the 2x2 case. However, the code seems clumsy. Is there any better way to do it?

set lang aql;
CREATE ARRAY InvertMe <val:double> [row=0:1,32,0, col=0:1,32,0];
INSERT INTO InvertMe '[[(1.75) (0.433013)] [(0.433013) (1.25)]]';

set lang afl;
load_library('dense_linear_algebra');

-- Singular value decomposition 
store(transpose(gesvd(InvertMe, 'U')), Ut);
store(project(apply(gesvd(InvertMe, 'S'),isigma,1/sigma),isigma),Splus);
store(transpose(gesvd(InvertMe, 'VT')), V);

-- Convert vector S+ to a diagonal matrix
store(project(apply(cross_join(V, Splus, V.col, Splus.i), nisigma, iif(col != i, 0, isigma)), nisigma), Splusmat);

-- Inverse matrix = V  *  S+  *  U'
store(gemm(gemm(V, Splusmat, build(V, 0)), Ut, build(V, 0)), Inverted);

-- Test the results
gemm(InvertMe, Inverted, build(InvertMe, 0));
gemm(Inverted, InvertMe, build(InvertMe, 0));

Kind regards.


#2

Hi. I showed your AFL snippet to one of our more high-powered math people, who had a question:

But… your code is basically doing what needs to be done; there is no more elegant shortcut.