This page needs Javascript
SINGULAR VALUE Decomposition SVD
Enter an NxM
matrix in the field 'Matrix A ' -
row by row, separating the elements with spaces /or tabs/. The resulted matrices
U
, S
and V
, such that A =
UxSxV
(P
is a Singular value
Decomposition ), will appear respectively in fields 'Matrix U ',
'Matrix S ' and 'Matrix V '. mailto:plarmuseau@pandora.be
Note: The javascript
does not work with INTEGERS ! .
Example with O-matrix /Mathlab
O>(u,s,v)=svd(x)
O>x
{
[ 1 , 2 ]
[ 2 , 3 ]
[ 3 , 4 ]
}
O>[u,s,v]=svd(x)
O>u
{
[ -0.338098 , 0.847952 , 0.408248 ]
[ -0.550649 , 0.173547 , -0.816497 ]
[ -0.763201 , -0.500857 , 0.408248 ]
}
O>s
{
[ 6.54676 , 0 ]
[ 0 , 0.374153 ]
[ 0 , 0 ]
}
O>v
{
[ -0.569595 , -0.821926 ]
[ -0.821926 , 0.569595 ]
}
O>[e,d]=eigen(x'x)
O>[e,d]=eigen(x'*x)
O>e
{
[ (0.569595,0) , (-0.821926,0) ]
[ (0.821926,0) , (0.569595,0) ]
}
O>d
{
[ (42.86,0) , (0,0) ]
[ (0,0) , (0.139991,0) ]
}
O>(1/sqrt(d(1,1)))*x*e
{
[ (0.338098,0) , (0.0484613,0) ]
[ (0.550649,0) , (0.00991838,0) ]
[ (0.7632,0) , (-0.0286245,0) ]
}
O>(1/sqrt(d(1,1)))*x*v
{
[ (-0.338098,0) , (0.0484613,0) ]
[ (-0.550649,0) , (0.00991837,0) ]
[ (-0.7632,0) , (-0.0286245,0) ]
}
O>(1/sqrt(d(2,2)))*x*d
{
[ (114.552,0) , (0.748306,0) ]
[ (229.104,0) , (1.12246,0) ]
[ (343.656,0) , (1.49661,0) ]
}
O>(1/sqrt(d(2,2)))*x*e
{
[ (5.91588,0) , (0.847952,0) ]
[ (9.635,0) , (0.173547,0) ]
[ (13.3541,0) , (-0.500858,0) ]
}
O>(1/sqrt(d(2,2)))*x*v
{
[ (-5.91588,0) , (0.847952,0) ]
[ (-9.635,0) , (0.173547,0) ]
[ (-13.3541,0) , (-0.500858,0) ]
}
O>x*v
{
[ -2.21345 , 0.317264 ]
[ -3.60497 , 0.0649331 ]
[ -4.99649 , -0.187398 ]
}
O>v
{
[ -0.402663 , 0.915348 ]
[ -0.915348 , -0.402663 ]
}
O>u
{
[ -0.32311 , 0.853776 , 0.408248 ]
[ -0.547507 , 0.18322 , -0.816497 ]
[ -0.771904 , -0.487337 , 0.408248 ]
}
O>1/s'
{
[ 0.24515 , 1.#INF , 1.#INF ]
[ 1.#INF , 1.6653 , 1.#INF ]
}
O>si
{
[ 0.24515 , 0 , 0 ]
[ 0 , 1.6653 , 0 ]
}
O>v*si*u'
{
[ 1.33333 , 0.333333 , -0.666665 ]
[ -0.499999 , 4.36201e-007 , 0.5 ]
}
O>pinv(x)
{
[ 1.33333 , 0.333333 , -0.666667 ]
[ -0.5 , -6.24924e-008 , 0.5 ]
}
O>svd(x)
{
4.07914
0.600491
}
External links and good applications
This is actually the most forgotten type of matrix math, although
certainly will prove to be the next decade the most important code wherewith
one will develop a myriad of applications
Weather
Forecast El Nino
Noise
filtering
Image
compression
Latent
Semantic Indexing: say creating Google type indexes Beside this
application, SVD could be used to OCR an text, to translate with a trained
matrix, to store knowledge, to create a knowledge matrix, it shaped the
internet indexing already, but a scaledown should change your computer or
database indexing, but think of homeselection, computer buying, book
suggestion, idee suggestor, optimal partner selection and so on
Source code is inspired from NUMERICAL RECIPES
IN C: THE ART OF SCIENTIFIC COMPUTING SVDCMP source code Integrated in
the javascript Check the matrix calculator
C++ source code
fp_t pythag(fp_t a, fp_t b)
{
float absa,absb;
absa=fabs(a);
absb=fabs(b);
if(absa>absb)
return absa*sqrt(1.0+sqr(absb/absa));
else
return(absb==0.0)?0.0:absb*sqrt(1.0+sqr(absa/absb));
}
void svdcmp(fp_t **a, int m, int n, fp_t *w, fp_t **v)
{
fp_t g,scale,anorm;
fp_t *rv1=new fp_t [n]-1;
g=scale=anorm=0.0;
for(int i=1;i<=n;++i){
int l=i+1;
rv1[i]=scale*g;
g=scale=0.0;
if(i<=m){
for(int k=i;k<=m;++k) scale+=fabs(a[k][i]);
if(scale){
fp_t s=0.0;
for(int k=i;k<=m;++k){
a[k][i]/=scale;
s+=a[k][i]*a[k][i];
}
fp_t f=a[i][i];
g=-sign(sqrt(s),f);
fp_t h=f*g-s;
a[i][i]=f-g;
for(int j=l;j<=n;++j){
fp_t sum=0.0;
for(int k=i;k<=m;++k) sum+=a[k][i]*a[k][j];
fp_t fct=sum/h;
for(int k=i;k<=m;++k) a[k][j]+=fct*a[k][i];
}
for(int k=i;k<=m;++k) a[k][i]*=scale;
}
}
w[i]=scale*g;
g=scale=0.0;
if((i<=m)&&(i!=n)){
for(int k=l;k<=n;++k) scale+=fabs(a[i][k]);
if(scale){
fp_t s=0.0;
for(int k=l;k<=n;++k){
a[i][k]/=scale;
s+=a[i][k]*a[i][k];
}
fp_t f=a[i][l];
g=-sign(sqrt(s),f);
fp_t h=f*g-s;
a[i][l]=f-g;
for(int k=l;k<=n;++k) rv1[k]=a[i][k]/h;
for(int j=l;j<=m;++j){
fp_t sum=0.0;
for(int k=l;k<=n;++k) sum+=a[j][k]*a[i][k];
for(int k=l;k<=n;++k) a[j][k]+=sum*rv1[k];
}
for(int k=l;k<=n;++k) a[i][k]*=scale;
}
}
anorm=max(anorm,(fabs(w[i])+fabs(rv1[i])));
}
{
fp_t f;
for(int i=n,l;i>=1;--i){
if(i<n){ // this makes f and l not dependent
if(f){
for(int j=l;j<=n;++j) v[j][i]=(a[i][j]/a[i][l])/f;
for(int j=l;j<=n;++j){
fp_t sum=0.0;
for(int k=l;k<=n;++k) sum+=a[i][k]*v[k][j];
for(int k=l;k<=n;++k) v[k][j]+=sum*v[k][i];
}
}
for(int j=l;j<=n;++j) v[i][j]=v[j][i]=0.0;
}
v[i][i]=1.0;
f=rv1[i];
l=i;
}
}
for(int i=min(m,n);i>=1;--i){
int l=i+1;
g=w[i];
for(int j=l;j<=n;++j) a[i][j]=0.0;
if(g){
g=1.0/g;
for(int j=l;j<=n;++j){
fp_t sum=0.0;
for(int k=l;k<=m;++k) sum+=a[k][i]*a[k][j];
fp_t f=(sum/a[i][i])*g;
for(int k=i;k<=m;++k) a[k][j]+=f*a[k][i];
}
for(int j=i;j<=m;++j) a[j][i]*=g;
}else for(int j=i;j<=m;++j) a[j][i]=0.0;
++a[i][i];
}
for(int k=n;k>=1;--k){
for(int its=1;its<=30;++its){
int flag=1,nm,l;
for(l=k;l>=1;--l){
nm=l-1;
if((fp_t)(fabs(rv1[l])+anorm)==anorm){
flag=0;
break;
}
if((fp_t)(fabs(w[nm])+anorm)==anorm) break;
}
if(flag){
fp_t c=0.0,s=1.0;
for(int i=l;i<=k;++i){
fp_t f=s*rv1[i];
rv1[i]=c*rv1[i];
if((fp_t)(fabs(f)+anorm)==anorm) break;
g=w[i];
fp_t h=pythag(f,g);
w[i]=h;
h=1.0/h;
c=g*h;
s=-f*h;
for(int j=1;j<=m;++j){
fp_t y=a[j][nm];
fp_t z=a[j][i];
a[j][nm]=y*c+z*s;
a[j][i]=z*c-y*s;
}
}
}
fp_t z=w[k];
if(l==k){
if(z<0.0){
w[k]=-z;
for(int j=1;j<=n;++j) v[j][k]=-v[j][k];
}
break;
}
if(its==30) exit(fprintf(stderr,"no convergence in 30 svdcmp iterations"));
fp_t x=w[l];
fp_t y=w[nm=k-1];
g=rv1[nm];
fp_t h=rv1[k];
fp_t f=((y-z)*(y+z)+(g-h)*(g+h))/(2.0*h*y);
g=pythag(f,1.0);
f=((x-z)*(x+z)+h*((y/(f+sign(g,f)))-h))/x;
fp_t c=1.0,s=1.0;
for(int j=l;j<=nm;++j){
int i=j+1;
g=rv1[i];
y=w[i];
h=s*g;
g*=c;
z=pythag(f,h);
rv1[j]=z;
c=f/z;
s=h/z;
f=x*c+g*s;
g=g*c-x*s;
h=y*s;
y*=c;
for(int jj=1;jj<=n;++jj){
x=v[jj][j];
z=v[jj][i];
v[jj][j]=x*c+z*s;
v[jj][i]=z*c-x*s;
}
z=pythag(f,h);
w[j]=z;
if(z){
z=1.0/z;
c=f*z;
s=h*z;
}
f=c*g+s*y;
x=c*y-s*g;
for(int jj=1;jj<=m;++jj){
y=a[jj][j];
z=a[jj][i];
a[jj][j]=y*c+z*s;
a[jj][i]=z*c-y*s;
}
}
rv1[l]=0.0;
rv1[k]=f;
w[k]=x;
}
}
delete (rv1+1);
}