Solved

Complex SVD algorithm

Posted on 1997-03-09
6
2,222 Views
Last Modified: 2010-04-16
I'm implementing the SVD with complex class elements.
So far I've modified the svdcmp() from the book "Numerical Recipes for C" by Press et al. for complex elements as outlined in the LINPACK user manual.  The bidiagonal form now has real elements (w[i] and rv1[i], and its appropriate Householder transforms in U and V matrices).  But the accumulation of the real rotations (Diagonalization part) in the _complex_ matrices U and V does not give the right results.  "Numerical ..." does not give explanation of its svdcmp(), so I wonder if anyone can explain how the accummulation of the rotations are done in the U and V matrices.
Thanks a huge bunch.
PK

Looks like I've to upload the source.  Don't have the source right now; will upload it soon.  By the way, this question might be harder than its current worth.  
Thanks again.
PK


Here's the _incomplete_ source:
Note before reading:
1. I've modified it half-way only.
2. To use complex arithmetic in the Householder part is
   beyond my capability.
3. Anyway, I've just implemented the SVD by breaking the
   complex matrices down to real matrices and process them
   just like normal.  Only thing is that the    number of operations have increased by two.  But I'm    working with small matrices, so it's okay.
4. Thanks for the help guys, how do I clear this question
   anyway?
   

#include <stdio.h>
#include <stddef.h>
#include <stdlib.h>
#include <conio.h>
#define NR_END 1
#include <math.h>
#include <complex.h>    // I used the complex C++ class
                        // in BC 5.0
void nrerror(char error_text[])
{
   fprintf(stderr,"Numerical Recipes run-time error...\n");
   fprintf(stderr,"%s\n",error_text);
   fprintf(stderr,"...now exiting to system...\n");
   exit(1);
}

complex *cpvector(long nl, long nh)
{
   complex *v;
   long n=nh-nl+1;
   v=new complex[n+NR_END];
   return v;
}

void free_cpvector(complex *v, long nl, long nh)
{
   delete[] v;
}

double *dvector(long nl, long nh)
{
   double *v;
   long n=nh-nl+1;
   v=new double[n+NR_END];
   return v;
}

void free_dvector(double *v, long nl, long nh)
{
   delete[] v;
}

complex **cpmatrix(long nrl, long nrh, long ncl, long nch)
{
   long i, j, nrow=nrh-nrl+1, ncol=nch-ncl+1;
   complex **m;

   /* allocate pointers to rows */
   m = new complex*[nrow+NR_END];
   for (i = 0; i < nrow; i++)
      m[i] = new complex[ncol+NR_END];
   return m;
}

void free_cpmatrix(complex **m, long nrl, long nrh, long ncl, long nch)
{
   int i, nrow = nrh-nrl;

   for (int i = 0; i < nrow+1;  i++)
      delete[] m[i];
   delete[] m;
}

void svdcmp(complex **a, int m, int n, complex w[], complex **v)
{
   double pythag(complex a, complex b);
   int flag, i, its, j, jj, k, l, nm;
   double anorm;
   complex c, f, g, h, s, scale, x, y, z, *rv1;
   double  C, F, G, H, S, X, Y, Z, *W, *RV1;
   rv1=cpvector(1,n);
   RV1=dvector(1,n);
   W=dvector(1,n);
   B=cpmatrix(1,m,1,n);
   g=scale=complex(0.0, 0.0);
   anorm=0.0;

   for (i=1;i<=n;i++) {
      l=i+1;
      rv1[i]=scale*g;
      g=s=scale=complex(0.0, 0.0);
      if (i <= m) {
         for(k=i;k<=m;k++) scale += complex(abs(a[k][i]),0);
            if (scale != complex(0.0,0.0)) {
               for (k=i;k<=m;k++) {
                  a[k][i] /= scale;
                  s += a[k][i]*a[k][i];
               }
               f=a[i][i];
               //g=-SIGN(sqrt(s),f);
               if (f != complex(0.0,0.0))
                  g=-sqrt(s)*f/abs(f);
               else g=-sqrt(s);                                                
               h=f*g-s;
               a[i][i]=f-g;
               for(j=l;j<=n;j++) {
                  for(s=complex(0.0,0.0),k=i;k<=m;k++) s += a[k][i]*a[k][j];
               f=s/h;
               for(k=i;k<=m;k++) a[k][j] += f*a[k][i];
            }
            for (k=i;k<=m;k++) a[k][i] *= scale;
         }
      }
      w[i]=scale*g;
      g=s=scale=complex(0.0,0.0);
      if (i<=m && i!=n) {
        for(k=l;k<=n;k++) scale += complex(abs(a[i][k]),0.0);
         if (scale != complex(0.0,0.0)) {
            for (k=l;k<=n;k++) {
               a[i][k] /= scale;
               s += a[i][k]*a[i][k];
            }
            f=a[i][l];
            //g=-SIGN(sqrt(s),f);
            if (f != complex(0.0,0.0))
               g=-sqrt(s)*f/abs(f);
            else g=-sqrt(s);
            h=f*g-s;
            a[i][l]=f-g;
            for (k=l;k<=n;k++) rv1[k]=a[i][k]/h;
            for (j=l;j<=m;j++) {
               for(s=complex(0.0,0.0),k=l;k<=n;k++)
                  s +=  a[j][k]*a]i][k];
               for (k=l;k<=n;k++) a[j][k] += s*rv1[k];
            }
            for (k=l;k<=n;k++) a[i][k] *= scale;
         }
      }
      anorm=FMAX(anorm, (abs(w[i])+abs(rv1[i])));
   }

   for(i=n;i>=1;i--){  // Accumulation of RH transformations
      if (i<n) {
         if (g != complex(0.0,0.0)) {
            for(j=l;j<=n;j++)
              v[j][i]=(a[i][j]/a[i][l])/g; //Double division
               for(j=l;j<=n;j++) {
               for(s=complex(0.0,0.0),k=l;k<=n;k++)
                  s += a[i][k]*v[k][j];
               for(k=l;k<=n;k++) v[k][j] += s*v[k][i];
            }
         }
         for(j=l;j<=n;j++) v[i][j]=v[j][i]=complex(0.0,0.0);
      }
      v[i][i]=complex(1.0,0.0);
      g=rv1[i];
      l=i;
   }

   for(i=IMIN(m,n);i>=1;i--){ //Accumulation of LH                                   transformations
      l=i+1;
      g=w[i];
      for(j=l;j<=n;j++) a[i][j]=complex(0.0,0.0);
      if (g != complex(0.0,0.0)) {
         g=1.0/g;
         for(j=l;j<=n;j++) {
            for(s=complex(0.0,0.0),k=l;k<=m;k++)
               s += a[k][i]*a[k][j];
            f=(s/a[i][i])*g;
               for(k=i;k<=m;k++) a[k][j] += f*a[k][i];
         }
         for(j=i;j<=m;j++) a[j][i] *= g;
      } else for(j=i;j<=m;j++) a[j][i]=complex(0.0,0.0);
      a[i][i]=1.0+a[i][i];
   }

   c=complex(1.0,0.0);
   for(i=1;i<n;i++) {
      if ((w[i] != complex(0.0,0.0))&&(rv1[i+1] !=                    complex(0.0,0.0))) {
         for(k=1;k<=m;k++) {
            a[k][i] *= (w[i]/abs(w[i]))*c;
            v[k][i+1] *=                (rv1[i+1]/abs(rv1[i+1]))*(abs(w[i])/w[i])/c;
         }
         c=(rv1[i+1]/abs(rv1[i+1]))*(abs(w[i])/w[i])/c;
      }
      else if ((w[i] == complex(0.0,0.0))&&(rv1[i+1] !=                         complex(0.0,0.0))) {
         for(k=1;k<=m;k++)
            v[k][i+1] *= (rv1[i+1]/abs(rv1[i+1]));
         c=(rv1[i+1]/abs(rv1[i+1]));
      }
      else if ((w[i] != complex(0.0,0.0))&&(rv1[i+1] ==                         complex(0.0,0.0))) {
         for(k=1;k<=m;k++) a[k][i] *= (w[i]/abs(w[i]))*c;
         c=complex(1.0,0.0);
      }
      else c=complex(1.0,0.0);
   }
   for(k=1;k<=m;k++) a[k][n] *= (w[n]/abs(w[n]))*c;

   for(i=1;i<=n;i++) w[i] = complex(W[i]=abs(w[i]),0.0);
   for(i=1;i<=n;i++)
      rv1[i] = complex(RV1[i]=abs(rv1[i]),0.0);

   printf("Matrix LH\n");
   for (k=1;k<=m;k++) {
      for (l=1;l<=n;l++)
            printf("%12.6f %12.6f", real(a[k][l]),                                 imag(a[k][l]));
      printf("\n");
   }

   printf("Matrix RH\n");
   for (k=1;k<=m;k++) {
      for (l=1;l<=n;l++)
            printf("%12.6f %12.6f", real(v[k][l]),                                 imag(v[k][l]));
      printf("\n");
   }

   for(k=n;k>=1;k--){  // Diagonalization of bidiagonal form
      for(its=1;its<=30;its++){ // Loop over singular values
         flag=1;
         for(l=k;l>=1;l--) {  // Test for splitting
            nm=l-1;           // Note rv1[1] is always zero
            if ((double)(abs(RV1[l])+anorm) == anorm) {
               flag=0;
               break;
            }
            if ((double)(abs(W[nm])+anorm) == anorm) break;
         }
         if (flag) {
            C=0.0;//c=complex(0.0,0.0);
            S=0.0;//s=complex(1.0,0.0);
            for(i=l;i<=k;i++) {
            F=S*RV1[i];
               RV1[i]=C*RV1[i];
               if ((double)(abs(F)+anorm) == anorm) break;
               G=W[i];
               H=pythag(F,G);
               W[i]=H;
               H=1.0/H;
               C=G*H;
               S=-F*H;
               for(j=1;j<=m;j++) {
                  y=a[j][nm]);
                  z=a[j][i]);
                  a[j][nm]=y*C+z*S;
                  a[j][i]=z*C-y*S;
               }
            }
         }
         Z=W[k];
         if (l==k) {     // Convergence
            if (Z<0.0) { // Singular value is made                                   non-negative
               W[k] = -Z;
               for(j=1;j<=n;j++) v[j][k] = -v[j][k];
            }
            break;
         }
         if (its == 30) nrerror("no convergence in 30 svdcmp                                      iterations");
         X=W[l];       // Shift from bottom 2 by 2 minor
         nm=k-1;
         Y=W[nm];
         G=RV1[nm];
         H=RV1[k];
         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;
         C=S=1.0;           // Next QR transformation
         for (j=l;j<=nm;j++) {
            i=j+1;
            G=RV1[i];
            Y=W[i];
            H=S*G;
            G=C*G;
            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(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;   // Rotation can be arbitrary if z=0
            if (Z != 0.0) {
               Z=1.0/Z;
               C=F*Z;
               S=H*Z;
            }
            F=C*G+S*Y;
            X=C*Y-S*G;
            for(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;
      }
   }
   for(i=1;i<=n;i++) w[i] = complex(W[i],0.0);
   for(i=1;i<=n;i++) rv1[i] = complex(RV1[i],0.0);
   free_dvector(W,1,n);
   free_dvector(RV1,1,n);
   free_cpvector(rv1,1,n);
}

double pythag(complex a, complex b)
{
   double absa, absb;
   absa=abs(a);
   absb=abs(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 svbksb(complex **u, complex w[], complex **v, int m, int n, complex b[], complex x[])
{
   int jj,j,i;
   complex s, *tmp;

   tmp=cpvector(1,n);
   for (j=1;j<=n;j++) {
      s=0.0;
      if (w[j] != complex(0.0,0.0)) {
         for (i=1;i<=m;i++) s += u[i][j]*b[i];
         s /= w[j];
      }
      tmp[j]=s;
   }
   for (j=1;j<=n;j++) {
      s=complex(0.0,0.0);
      for (jj=1;jj<=n;jj++) s += v[j][jj]*tmp[jj];
      x[j]=s;
   }
   free_cpvector(tmp,1,n);
}

void main(void)
{
   complex **A, **B, *W, **U, **V, *X, *C;
   int j,k,l,m=3,n=3, sm=2, sn=3;;
   double wmax,wmin;

   W=cpvector(1,5);
   X=cpvector(1,5);
   C=cpvector(1,5);
   A=cpmatrix(1,5,1,5);
   B=cpmatrix(1,5,1,5);
   U=cpmatrix(1,5,1,5);
   V=cpmatrix(1,5,1,5);
   A[1][1] = complex(1.0 , 0.0);
   A[1][2] = complex(2.0 , 0.0);
   A[1][3] = complex(3.0 , 0.0);
   A[2][1] = complex(2.0 , 0.0);
   A[2][2] = complex(2.0 , 0.0);
   A[2][3] = complex(3.0 , 0.0);
   A[3][1] = complex(3.0 , 0.0);
   A[3][2] = complex(3.0 , 0.0);
   A[3][3] = complex(3.0 , 0.0);
   B[1][1] = complex(1.0 , 0.0);
   B[2][1] = complex(1.0 , 0.0);
   B[3][1] = complex(1.0 , 0.0);
   B[1][2] = complex(1.0 , 0.0);
   B[2][2] = complex(2.0 , 0.0);
   B[3][2] = complex(3.0 , 0.0);
   for(k=1;k<=m;k++)
      for(l=1;l<=n;l++) U[k][l]=A[k][l];
   if (n>m) {
      for (k=m+1;k<=n;k++) {
         for (l=1;l<=n;l++) {
            A[k][l] = complex(0.0,0.0);
            U[k][l] = complex(0.0,0.0);
         }
      }
      m=n;
   }
   svdcmp(U,m,n,W,V);
   printf("Decomposition Matrices:\n");
   printf("Matrix A\n");
   for (k=1;k<=m;k++) {
      for (l=1;l<=n;l++)
         printf("%12.6f %12.6f", real(A[k][l]),                                  imag(A[k][l]));
      printf("\n");
   }
   printf("Matrix U\n");
   for (k=1;k<=m;k++) {
      for (l=1;l<=n;l++)
         printf("%12.6f %12.6f", real(U[k][l]),                                  imag(U[k][l]));
      printf("\n");
   }
   printf("Diagonal of matrix w\n");
   for (k=1;k<=n;k++)
      printf("%12.6f %12.6f\n", real(W[k]), imag(W[k]));
   printf("Matrix v-transpose\n");
   for (k=1;k<=n;k++) {
      for (l=1;l<=n; l++)
         printf("%12.6f %12.6f  ", real(V[l][k]),                                    imag(V[l][k]));
      printf("\n");
   }
   
   wmax=0.0;
   for (k=1;k<=sn;k++)
      if (abs(W[k]) > wmax) wmax=abs(W[k]);
   wmin=wmax*(1.0e-6);
   for (k=1;k<=sn;k++)
      if (abs(W[k]) < wmin) W[k]=complex(0.0,0.0);

   for (l=1;l<=sm;l++) {
      printf("\nVector number %2d\n",l);
      for (k=1;k<=sn;k++) C[k]=B[k][l];
      svbksb(U,W,V,sn,sn,C,X);
      printf(" solution vector is:\n");
      for (k=1;k<=sn;k++) printf("%12.6f %12.6f",real(X[k]),                                                 imag(X[k]));
      printf("\n original right-hand side vector:\n");
      for (k=1;k<=sn;k++) printf("%12.6f %12.6f",real(C[k]),                                                 imag(C[k]));
      printf("\n (matrix)*(sol'n vector):\n");
      for (k=1;k<=sn;k++) {
             C[k]=0.0;
         for (j=1;j<=sn;j++)
            C[k] += A[k][j]*X[j];
      }
      for (k=1;k<=sn;k++) printf("%12.6f %12.6f",real(C[k]),                                                 imag(C[k]));
      printf("\n");
   }

   free_cpmatrix(V,1,5,1,5);
   free_cpmatrix(U,1,5,1,5);
   free_cpmatrix(B,1,5,1,5);
   free_cpmatrix(A,1,5,1,5);
   free_cpvector(C,1,5);
   free_cpvector(X,1,5);
   free_cpvector(W,1,5);
}
0
Comment
Question by:PK030997
6 Comments
 
LVL 2

Expert Comment

by:lockhart
Comment Utility
You need to include your source or make it available, along with an explanation of the SVD algorithm and what you  are trying to accomplish by modifying it.

Mike
0
 

Author Comment

by:PK030997
Comment Utility
Edited text of question
0
 
LVL 3

Expert Comment

by:LucHoltkamp
Comment Utility
I assume you are talking about Householder QU-decomposition.
Why do you want to use Householder decomposition. If you want
to solve ordinary lineair equations you should use Gauss Pivoting, its much simpler and faster.
If however you still want to use it, for least squares or something you should buy a good book about mathematics or pay (much) more points for source. I have to warn you: the source will be C++
0
Free Trending Threat Insights Every Day

Enhance your security with threat intelligence from the web. Get trending threat insights on hackers, exploits, and suspicious IP addresses delivered to your inbox with our free Cyber Daily.

 

Author Comment

by:PK030997
Comment Utility
Edited text of question
0
 

Author Comment

by:PK030997
Comment Utility
Edited text of question
0
 
LVL 1

Accepted Solution

by:
roy020697 earned 200 total points
Comment Utility
I believe the only way to clear a question is to receive an answer (like this one) and give it a passing grade. If points are valuable to you then you can adjust the points to something small then edit your question to be easily answered. I guess your complex SVD algorithm question has been answered?  :-)
0

Featured Post

Threat Intelligence Starter Resources

Integrating threat intelligence can be challenging, and not all companies are ready. These resources can help you build awareness and prepare for defense.

Join & Write a Comment

Have you thought about creating an iPhone application (app), but didn't even know where to get started? Here's how: ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ Important pre-programming comments: I’ve never tri…
Summary: This tutorial covers some basics of pointer, pointer arithmetic and function pointer. What is a pointer: A pointer is a variable which holds an address. This address might be address of another variable/address of devices/address of fu…
The goal of this video is to provide viewers with basic examples to understand how to use strings and some functions related to them in the C programming language.
The goal of this video is to provide viewers with basic examples to understand how to create, access, and change arrays in the C programming language.

762 members asked questions and received personalized solutions in the past 7 days.

Join the community of 500,000 technology professionals and ask your questions.

Join & Ask a Question

Need Help in Real-Time?

Connect with top rated Experts

12 Experts available now in Live!

Get 1:1 Help Now