I found a better solution and implemented it.
The solution is based on the following. Compute the singular value
decomposition of the matrix, let this be U*S*V^T (U and V are
orthogonal, S is diagonal, ^T means transposition). Now x*U*U^T is
the projection of vector x to the row space of the matrix.
I use the GSL library to compute the singular value decomposition.
One could use LAPACK instead, which would have the advantage that
LAPACK has svd implemented for complex matrices as well.
In my question I asked the solution to work even if the matrix is not
of complete rank. This might not be a reasonable request, for such a
case is very instable in the matrix. This solution nevertheless
attempts to honor this request using tolerant comparisions of J.
Let's see the code -- C wrapper first, then J code.
$ cat svdwrap.c
#include <gsl/gsl_vector_double.h>
#include <gsl/gsl_matrix_double.h>
#include <gsl/gsl_linalg.h>
/*
wrap_svd_decompose(m, n, A, S, V)
Performs SVD decomposition on matrix of doubles.
In condition: m <= n.
A is input matrix, size m x n,
U is output matrix, same place as A, size m x n,
S is output vector, size n x n,
V is output matrix, size n x n.
Out condition: U and V are orthogonal matrices (ie. U*U^T=I),
S is all nonnegative and sorted non-increasing,
A = U*diag(S)*V^T
*/
int
wrap_svd_decompose(
int zm, int zn, double *ac, double *sc, double *vc
) {
if (zn < 0 || zm < zn)
return 1;
gsl_matrix_view ai = gsl_matrix_view_array(ac, zm, zn);
gsl_vector_view si = gsl_vector_view_array(sc, zn);
gsl_matrix_view vi = gsl_matrix_view_array(vc, zn, zn);
gsl_vector *wv = gsl_vector_alloc(zn);
int err = gsl_linalg_SV_decomp(&ai.matrix, &vi.matrix, &si.vector, wv);
gsl_vector_free(wv);
return err;
}
$ gcc -lgsl -lgslcblas -Wall -O -fPIC -shared -o svdwrap.so svdwrap.c
$ cat svd.ijs
NB. singular value decomposition
svd=: 3 :0"2
tr=. 0
if. 2>#$y do. 13!:8]3 end.
if. tr=. </$y do. y=. |: y end.
'mz nz'=: $y
p=. mz;nz;y;(i.nz);(i.,~nz)
'err mr nr ut si vt'=. './svdwrap.so wrap_svd_decompose i i i *d *d *d' 15!:0 p
if. err do. 13!:8]10 end.
if. tr do. vt;si;ut else. ut;si;vt end.
)
svd_tol=: 3 :0"2
'ut si vt'=.svd y
mb=. ({.<{.+])si
(mb#"1 ut);(mb#si);(mb#"1 vt)
)
NB. orthogonal projection
oproj=: 4 :0
ut=. >{.svd_tol |:y
(x +/ .* ut) +/ .* |: ut
)
$ j701console # j701_64 on amd64-linux
0!:0<'svd.ijs'
NB. test svd
NB.a=: 4 _9 _4 ,: 4 _3 _2
NB.a=: 0 1 2,:0 1 2
NB.a=: 10-?4 10$21
a=: 10-?10 4$21
'ut si vt'=: svd_tol a
NB.echo a;ut;si;vt
NB. following four numbers must be small
>./|,(=i.#{.ut)- ut+/ .*~|:ut
6.66134e_16
>./|,(=i.#{.vt)- vt+/ .*~|:vt
6.66134e_16
-.*./2>:/\si
0
>./|,a- ut+/ .*si*|:vt
7.10543e_15
NB. test orthogonal projection
ve=: _1 _7 12
sp=: 4 _9 _4 ,: 4 _3 _2
ve oproj sp
_4 _3 _4.71845e_16
sp=: 10-?4 10 $ 21
ve=: 10-?10 $ 21
NB. following three should be equal
ve oproj sp
1.08054 _6.34628 1.53816 _1.43236 0.630816 _4.42371 _1.54375 _1.3119
_4.15668 1.93827
ve (] +/ .* %.) |: sp
1.08054 _6.34628 1.53816 _1.43236 0.630816 _4.42371 _1.54375 _1.3119
_4.15668 1.93827
ve ([-[:{.[:((],~[-[:+/+/ .*~*(%+/&:*:"1)@]),:^:(2>#...@$))/,) sp
1.08054 _6.34628 1.53816 _1.43236 0.630816 _4.42371 _1.54375 _1.3119
_4.15668 1.93827
exit 0
$
--
Ambrus
----------------------------------------------------------------------
For information about J forums see http://www.jsoftware.com/forums.htm
