search

optimatika / ojAlgo

oj! Algorithms
http://ojalgo.org
MIT License
459 stars 208 forks source link

Singular Value Decomposition computing wrong - Or am I using wrong options? #253

Closed DanielMartensson closed 4 years ago

DanielMartensson commented 4 years ago

Here is a minimal example of the SVD code

public class Svd {

    static Logger logger = LoggerFactory.getLogger(Svd.class);

    static public void svd(double A[][], double U[][], double S[], int rows, int columns) {

        // Create data
        Primitive64Matrix data = Primitive64Matrix.FACTORY.rows(A);

        // Print data
        for(int i = 0; i < rows; i++) {
            for(int j = 0; j < columns; j++) {
                System.out.print("\t" + data.get(i, j));
            }
            System.out.println("");
        }

        // Perform SVD
        SingularValue<Double> svd = SingularValue.PRIMITIVE.make();
        svd.decompose(data);
        MatrixStore<Double> svdU = svd.getU();
        MatrixStore<Double> svdS = svd.getD();
        MatrixStore<Double> svdV = svd.getV();
        logger.info("Done with Singular Value Decomposition");

        // Copy over to U, S
        for(int i = 0; i < rows; i++) {
            for(int j = 0; j < columns; j++) {
                U[i][j] = svdU.get(i, j);
                System.out.print("\t" + U[i][j]);
            }
            System.out.println("");
        }
        for(int i = 0; i < columns; i++) {
            S[i] = svdS.get(i, i);
            System.out.println(S[i]);
        }   
    }
}

The output looks like this:

    -3.0    -2.0    -1.0    0.0 1.0 2.0 3.0
    -3.0    -2.0    -1.0    0.0 1.0 2.0 3.0
    -3.0    -2.0    -1.0    0.0 1.0 2.0 3.0
    -3.0    -2.0    -1.0    0.0 1.0 2.0 3.0
    -3.0    -2.0    -1.0    0.0 1.0 2.0 3.0
    -3.0    -2.0    -1.0    0.0 1.0 2.0 3.0
    -3.0    -2.0    -1.0    0.0 1.0 2.0 3.0
    -3.0    -2.0    -1.0    0.0 1.0 2.0 3.0
    -3.0    -2.0    -1.0    0.0 1.0 2.0 3.0
    -3.0    -2.0    -1.0    0.0 1.0 2.0 3.0
[main] INFO se.danielmartensson.fisherfaces.matlab.ojalgo.Svd - Done with Singular Value Decomposition
    -0.316227766016838  0.9486832980505135  -1.1102230246251565E-16 -1.3877787807814457E-17 -4.1633363423443376E-17 -1.3877787807814457E-17 2.7755575615628914E-17
    -0.31622776601683794    -0.10540925533894596    0.9428090415820636  -3.334181172154923E-18  -1.0002543516464769E-17 -1.721196897996938E-17  -2.1087213271319068E-17
    -0.31622776601683794    -0.10540925533894602    -0.11785113019775786    0.9354143466934853  1.7040869798512316E-16  -6.803628124108537E-18  -2.7053155959738237E-19
    -0.31622776601683794    -0.10540925533894602    -0.11785113019775792    -0.13363062095621228    0.9258200997725512  2.1524097680092275E-16  2.7485044056031528E-17
    -0.31622776601683794    -0.10540925533894602    -0.11785113019775792    -0.13363062095621217    -0.15430334996209183    0.9128709291752767  2.0163162690960059E-16
    -0.31622776601683794    -0.10540925533894602    -0.11785113019775792    -0.13363062095621217    -0.15430334996209188    -0.18257418583505536    0.8944271909999159
    -0.31622776601683794    -0.10540925533894602    -0.11785113019775792    -0.13363062095621217    -0.15430334996209188    -0.18257418583505536    -0.223606797749979
    -0.31622776601683794    -0.10540925533894602    -0.11785113019775792    -0.13363062095621217    -0.15430334996209188    -0.18257418583505536    -0.223606797749979
    -0.31622776601683794    -0.10540925533894602    -0.11785113019775792    -0.13363062095621217    -0.15430334996209188    -0.18257418583505536    -0.223606797749979
    -0.31622776601683794    -0.10540925533894602    -0.11785113019775792    -0.13363062095621217    -0.15430334996209188    -0.18257418583505536    -0.223606797749979
16.733200530681515
0.0
0.0
0.0
0.0
0.0
0.0

But according to Octave, the output looks like this:

[U, D, V] = svd(X, 'econ');
X
U
D
------------------------------------
X =

  -3  -2  -1   0   1   2   3
  -3  -2  -1   0   1   2   3
  -3  -2  -1   0   1   2   3
  -3  -2  -1   0   1   2   3
  -3  -2  -1   0   1   2   3
  -3  -2  -1   0   1   2   3
  -3  -2  -1   0   1   2   3
  -3  -2  -1   0   1   2   3
  -3  -2  -1   0   1   2   3
  -3  -2  -1   0   1   2   3

U =

  -0.31623   0.94868  -0.00000   0.00000  -0.00000   0.00000   0.00000
  -0.31623  -0.10541  -0.33333  -0.33333  -0.33333  -0.33333  -0.33333
  -0.31623  -0.10541   0.91667  -0.08333  -0.08333  -0.08333  -0.08333
  -0.31623  -0.10541  -0.08333   0.91667  -0.08333  -0.08333  -0.08333
  -0.31623  -0.10541  -0.08333  -0.08333   0.91667  -0.08333  -0.08333
  -0.31623  -0.10541  -0.08333  -0.08333  -0.08333   0.91667  -0.08333
  -0.31623  -0.10541  -0.08333  -0.08333  -0.08333  -0.08333   0.91667
  -0.31623  -0.10541  -0.08333  -0.08333  -0.08333  -0.08333  -0.08333
  -0.31623  -0.10541  -0.08333  -0.08333  -0.08333  -0.08333  -0.08333
  -0.31623  -0.10541  -0.08333  -0.08333  -0.08333  -0.08333  -0.08333

D =

Diagonal Matrix

   16.73320          0          0          0          0          0          0
          0    0.00000          0          0          0          0          0
          0          0   -0.00000          0          0          0          0
          0          0          0    0.00000          0          0          0
          0          0          0          0    0.00000          0          0
          0          0          0          0          0    0.00000          0
          0          0          0          0          0          0    0.00000

As you can see. the U matrix is not correct compared with Octave. Why? Is it because I'm using PRIMITIVE? Don't know what it means.

I have also a question about eigendecomposition. Here is my example code:

public class Eig {

    static Logger logger = LoggerFactory.getLogger(Eig.class);

    // A*D = D*B*V - Find D and V
    static public void eig(double A[][], double B[][], double D[], double V[][], int rows) {

        // Create eigA and eigB
        Primitive64Matrix eigA = Primitive64Matrix.FACTORY.rows(A);
        Primitive64Matrix eigB = Primitive64Matrix.FACTORY.rows(B);

        // Find inverse of eigB - Or pseudo inverse if B is singular
        eigB = eigB.invert();

        // Multiply eigB*eigA
        Primitive64Matrix data = eigB.multiply(eigA);

        // Perform eigendecomposition
        Eigenvalue<Double> eig = Eigenvalue.PRIMITIVE.make(data,false);
        boolean success = eig.decompose(data);
        if(success == false)
            logger.error("Could not perform eigenvalue decomposition!");
        MatrixStore<Double> Deig = eig.getD();
        MatrixStore<Double> Veig = eig.getV();

        // Copy over to D, V
        for(int i = 0; i < rows; i++) {
            D[i] = Deig.get(i, i);
            for(int j = 0; j < rows; j++) {
                V[i][j] = Veig.get(i, j);
            }
        }   
    }
}

I call that code with this data and get the output:

// Create random data
        int rows = 5;
        double[][] A = new double[rows][rows];
        double[][] B = {{0.7868535918842593,    0.25206982118540255,    0.8502514736661777, 0.957057138020801,  0.057816696367313236},
                {0.13879101503702707,   0.7128202999301166, 0.1301621659119453, 0.31855375514373074,    0.5217887100001347},
                {0.32528138064626044,   0.375732472964022,  0.8722279148638662, 0.2623839802297663, 0.6554277853039092},
                {0.29653713361927037,   0.5156537331156116, 0.0769602762384829, 0.5727731028190134, 0.22446173846280937},
                {0.48182976689909873,   0.19256889992091686,    0.6541198727551285, 0.6268943810659438, 0.42040774186487684}};
        for(int i = 0; i < rows; i++) {
            for(int j = 0; j < rows; j++) {
                A[i][j] = i+j;
            }
        }
        System.out.println("Print A");
        for(int i = 0; i < rows; i++) {
            for(int j = 0; j < rows; j++) {
                System.out.print("\t" + A[i][j]);
            }
            System.out.println("");
        }
        System.out.println("Print B");
        for(int i = 0; i < rows; i++) {
            for(int j = 0; j < rows; j++) {
                System.out.print("\t" + B[i][j]);
            }
            System.out.println("");
        }

        // Perform eig
        double[][] V = new double[rows][rows];
        double[] D = new double[rows];
        Eig.eig(A, B, D, V, rows);

        System.out.println("Print eigenvalues");
        for(int i = 0; i < rows; i++) {
            System.out.println(D[i]);
        }   
        System.out.println("Print eigenvectors");
        for(int i = 0; i < rows; i++) {
            for(int j = 0; j < rows; j++) {
                System.out.print("\t" + V[i][j]);
            }
            System.out.println("");
        }
Print A
    0.0 1.0 2.0 3.0 4.0
    1.0 2.0 3.0 4.0 5.0
    2.0 3.0 4.0 5.0 6.0
    3.0 4.0 5.0 6.0 7.0
    4.0 5.0 6.0 7.0 8.0
Print B
    0.7868535918842593  0.25206982118540255 0.8502514736661777  0.957057138020801   0.057816696367313236
    0.13879101503702707 0.7128202999301166  0.1301621659119453  0.31855375514373074 0.5217887100001347
    0.32528138064626044 0.375732472964022   0.8722279148638662  0.2623839802297663  0.6554277853039092
    0.29653713361927037 0.5156537331156116  0.0769602762384829  0.5727731028190134  0.22446173846280937
    0.48182976689909873 0.19256889992091686 0.6541198727551285  0.6268943810659438  0.42040774186487684
Print eigenvalues
-72.24303468512798
19.68243960526661
-1.263846900438926E-13
2.1790180928485528E-14
-2.598102440286286E-14
Print eigenvectors
    -0.843534098386463  6.888766304572221   -0.4103808547114451 -0.15064677005171725    0.5300697841954638
    -0.021650343361678242   -0.7125750738079105 0.8421978260532967  0.5003740977706913  0.5722830367666923
    0.2626516887639604  -2.49185111743384   -0.18408075846905578    -0.7342555017208549 -1.4339375479081005
    0.45288497844928716 -4.075205813120237  -0.5169085423758242 0.5699757903365146  -0.96925315126561
    -0.11785912852251558    0.9110743367448255  0.26917232950303754 -0.18544761633463291    1.3008378782115613

But GNU Octave shows

V =

  -0.8435341  -0.8140592   0.5179533   0.4271032   0.0036528
  -0.0216503   0.0842064  -0.7341037  -0.0259704  -0.1009847
   0.2626517   0.2944670  -0.1408616  -0.7128079   0.5296932
   0.4528850   0.4815752   0.4122211  -0.2048857  -0.7710436
  -0.1178591  -0.1076635  -0.0552091   0.5165609   0.3386822

D =

Diagonal Matrix

  -7.2243e+01            0            0            0            0
            0   1.9682e+01            0            0            0
            0            0   1.2203e-12            0            0
            0            0            0   1.0295e-14            0
            0            0            0            0  -6.6871e-14

From the same data

A = [0.0 1.0    2.0 3.0 4.0
    1.0 2.0 3.0 4.0 5.0
    2.0 3.0 4.0 5.0 6.0
    3.0 4.0 5.0 6.0 7.0
    4.0 5.0 6.0 7.0 8.0];

B = [0.7868535918842593 0.25206982118540255 0.8502514736661777  0.957057138020801   0.057816696367313236
    0.13879101503702707 0.7128202999301166  0.1301621659119453  0.31855375514373074 0.5217887100001347
    0.32528138064626044 0.375732472964022   0.8722279148638662  0.2623839802297663  0.6554277853039092
    0.29653713361927037 0.5156537331156116  0.0769602762384829  0.5727731028190134  0.22446173846280937
    0.48182976689909873 0.19256889992091686 0.6541198727551285  0.6268943810659438  0.42040774186487684];

% A*D = D*B*V - Find V and D
[V, D] = eig(pinv(B)*A)

Notice that I'm using random matrices.

apete commented 4 years ago

The eigenvalues and the singular values are the same, aren't they? The eigenvectors and/or the singular vectors are not unique so you cannot require them to be the same as with some other tool.

DanielMartensson commented 4 years ago

The eigenvalues and the singular values are the same, aren't they? The eigenvectors and/or the singular vectors are not unique so you cannot require them to be the same as with some other tool.

But the data is the same?

apete commented 4 years ago

Which data? Don't know what you mean.

DanielMartensson commented 4 years ago

Which data. Don't know what you mean.

Let's take the eigenvalue example first. I assume that solving that problem, the SVD problem would be much easier to solve later.

Here is the matlab code:

A = [0.0 1.0    2.0 3.0 4.0
    1.0 2.0 3.0 4.0 5.0
    2.0 3.0 4.0 5.0 6.0
    3.0 4.0 5.0 6.0 7.0
    4.0 5.0 6.0 7.0 8.0];

B = [0.7868535918842593 0.25206982118540255 0.8502514736661777  0.957057138020801   0.057816696367313236
    0.13879101503702707 0.7128202999301166  0.1301621659119453  0.31855375514373074 0.5217887100001347
    0.32528138064626044 0.375732472964022   0.8722279148638662  0.2623839802297663  0.6554277853039092
    0.29653713361927037 0.5156537331156116  0.0769602762384829  0.5727731028190134  0.22446173846280937
    0.48182976689909873 0.19256889992091686 0.6541198727551285  0.6268943810659438  0.42040774186487684];

% A*D = D*B*V - Find V and D
[V, D] = eig(pinv(B)*A)

A and B are data.

They will output this data:

V =

  -0.8435341  -0.8140592   0.5179533   0.4271032   0.0036528
  -0.0216503   0.0842064  -0.7341037  -0.0259704  -0.1009847
   0.2626517   0.2944670  -0.1408616  -0.7128079   0.5296932
   0.4528850   0.4815752   0.4122211  -0.2048857  -0.7710436
  -0.1178591  -0.1076635  -0.0552091   0.5165609   0.3386822

D =

Diagonal Matrix

  -7.2243e+01            0            0            0            0
            0   1.9682e+01            0            0            0
            0            0   1.2203e-12            0            0
            0            0            0   1.0295e-14            0
            0            0            0            0  -6.6871e-14

Compare that with my code above, it will not give the same result.

apete commented 4 years ago

Are you referring to the fact that the V matrices are different? Eigenvectors are not unique which means the V matrix is not unique.

DanielMartensson commented 4 years ago

Are you referring to the fact that the V matrices are different? Eigenvectors are not unique which means the V matrix is not unique.

Yes. For eig - V is different. For svd - U is different.

apete commented 4 years ago

Again, eigenvectors are not uniquely defined. https://stackoverflow.com/questions/13041178/could-we-get-different-solutions-for-eigenvectors-from-a-matrix

DanielMartensson commented 4 years ago

Again, eigenvectors are not uniquely defined. https://stackoverflow.com/questions/13041178/could-we-get-different-solutions-for-eigenvectors-from-a-matrix

So you mean that if eigenvectors can differ depending on what software and computer I'm using?

apete commented 4 years ago

Depending on which software, yes. Eigenvalues are uniquely defined, but not vectors.

On 14 May 2020, at 16:23, Daniel Mårtensson notifications@github.com wrote:

Again, eigenvectors are not uniquely defined. https://stackoverflow.com/questions/13041178/could-we-get-different-solutions-for-eigenvectors-from-a-matrix

So you mean that if eigenvectors can differ depending on what software and computer I'm using?

— You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

DanielMartensson commented 4 years ago

Depending on which software, yes. Eigenvalues are uniquely defined, but not vectors. … On 14 May 2020, at 16:23, Daniel Mårtensson @.***> wrote: Again, eigenvectors are not uniquely defined. https://stackoverflow.com/questions/13041178/could-we-get-different-solutions-for-eigenvectors-from-a-matrix So you mean that if eigenvectors can differ depending on what software and computer I'm using? — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

Good to know. I assume that it was no issue after all?

Normaly, my MATLAB code want to solve: [V, D] = eig(A, B)

And I know that OjAlgo can solve generalized eigenvalue problem. But your examples shows that A and B must be symmetrical. Is that a requirement for generalized eigenvalue problem?

apete commented 4 years ago

It’s a requirement for what what ojAlgo can handle at least. Essentially ojAlgo can handle the cases described here: https://www.netlib.org/lapack/lug/node54.html

On 14 May 2020, at 16:31, Daniel Mårtensson notifications@github.com wrote:

Depending on which software, yes. Eigenvalues are uniquely defined, but not vectors. … On 14 May 2020, at 16:23, Daniel Mårtensson @.***> wrote: Again, eigenvectors are not uniquely defined. https://stackoverflow.com/questions/13041178/could-we-get-different-solutions-for-eigenvectors-from-a-matrix So you mean that if eigenvectors can differ depending on what software and computer I'm using? — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

Good to know. I assume that it was no issue after all?

Normaly, my MATLAB code want to solve: [V, D] = eig(A, B)

And I know that OjAlgo can solve generalized eigenvalue problem. But your examples shows that A and B must be symmetrical. Is that a requirement for generalized eigenvalue problem?

— You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

DanielMartensson commented 4 years ago

It’s a requirement for what what ojAlgo can handle at least. Essentially ojAlgo can handle the cases described here: https://www.netlib.org/lapack/lug/node54.html

I understand. Well, perhaps I need to use [V, D] = eig(pinv(B)*A)

instead of [V, D] = eig(A, B)

Because only B is symmetrical. I'm not talking about MATLAB's routines here. I'm only using MATLAB code to describe math easily.

Thank you for the information.

DanielMartensson commented 4 years ago

It’s a requirement for what what ojAlgo can handle at least. Essentially ojAlgo can handle the cases described here: https://www.netlib.org/lapack/lug/node54.html … On 14 May 2020, at 16:31, Daniel Mårtensson @.> wrote: Depending on which software, yes. Eigenvalues are uniquely defined, but not vectors. … On 14 May 2020, at 16:23, Daniel Mårtensson @.> wrote: Again, eigenvectors are not uniquely defined. https://stackoverflow.com/questions/13041178/could-we-get-different-solutions-for-eigenvectors-from-a-matrix So you mean that if eigenvectors can differ depending on what software and computer I'm using? — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe. Good to know. I assume that it was no issue after all? Normaly, my MATLAB code want to solve: [V, D] = eig(A, B) And I know that OjAlgo can solve generalized eigenvalue problem. But your examples shows that A and B must be symmetrical. Is that a requirement for generalized eigenvalue problem? — You are receiving this because you commented. Reply to this email directly, view it on GitHub, or unsubscribe.

En sak! Nu är denna tråd redan stängd så jag tar det på Svenska. Jag håller på skriva ett javaprogram som håller på med bildigenkänning med minimal data. Hittils har det visat sig goda resultat.

Jag skriver av ett 8-årigt MATLAB projekt som håller på med något som heter Fisherfaces. Det ska tydligen vara ett krångligt ämne. Men så länge jag skriver av MATLAB koden så blir det fint.

Men då finns det ett problem. MATLAB projektet använder sig utav egenvektorer från EIG och egenvektorer från SVD som multipliceras med varandra. Jag har ingen aning om hur matematiken går till då det inte finns beskrivet.

I MATLAB så blir det alltså bättre resulat (96%) men i mitt program så blir det 93.3% i validering. Nog är detta högt. Men jag misstänker att programmet är konstruerat för just MATLAB's SVD U och EIG V.

Finns det möjlighet inom OjAlgo att välja om så det blir "MATLAB-liknande" EIG V och SVD U? Typ någon val utav norm?

Det som skiljer mellan SVD U i MATLAB och SVD U i OjAlgo är att vissa värden är negativa, andra är positiva. Men själva nummret är det samma.

I EIG V så skiljer nummret och om det är positivt eller negativt.

Koden kommer jag lägga upp inom kort.