Search code examples
c++linear-algebraqr-codeprojectionorthogonal

Orthogonalization in QR Factorization outputting slightly innaccurate orthogonalized matrix


I am writing code for QR Factorization and for some reason my orthogonal method does not work as intended. Basically, my proj() method is outputting random projections. Here is the code:

apmatrix<double> proj(apmatrix<double> v, apmatrix<double> u)   
//Projection of u onto v
{
//proj(v,u) = [(u dot v)/(v dot v)]*v
    double a = mult(transpose(u,u),v)[0][0], b = mult(transpose(v,v),v)[0][0], c = (a/b);
    apmatrix<double>k;
    k.resize(v.numrows(),v.numcols());
    for(int i = 0; i<v.numrows(); i++)
    {
        for(int j = 0; j<v.numcols(); j++)
        {
            k[i][j]=v[i][j]*c;
        }
    }
    return k;
}

I tested the method by itself with manual matrix inputs, and it seems to work fine. Here is my orthogonal method:

apmatrix<double> orthogonal(apmatrix<double> A)   //Orthogonal
{
    /*
    n = (number of columns of A)-1
    x = columns of A
    v0 = x0
    v1 = x1 - proj(v0,x1)
    vn = xn - proj(v0,xn) - proj(v1,xn) - ... - proj(v(n-1),xn)
    V = {v1, v2, ..., vn} or [v0 v1 ... vn]
    */
    apmatrix<double> V, x, v;
    int n = A.numcols();
    V.resize(A.numrows(),n);
    x.resize(A.numrows(), 1);
    v.resize(A.numrows(),1);
    for(int i = 0; i<A.numrows(); i++)
    {
        x[i][0]=A[i][1];
        v[i][0]=A[i][0];
        V[i][0]=A[i][0];
    }
    for (int c = 1; c<n; c++)   //Iterates through each col of A as if each was its own matrix
    {
        apmatrix<double>vn,vc; //vn = Orthogonalized v (avoiding matrix overwriting of v); vc = previously orthogonalized v
            vn=x;
        vc.resize(v.numrows(), 1);
        for(int i=0; i<c; i++)   //Vn = an-(sigma(t=1, n-1, proj(vt, xn))
        {
            for(int k = 0; k<V.numrows(); k++)
                vc[k][0] = V[k][i]; //Sets vc to designated v matrix
            apmatrix<double>temp = proj(vc, x);
            for(int j = 0; j<A.numrows(); j++)
            {
                vn[j][0]-=temp[j][0]; //orthogonalize matrix
            }
        }
        for(int k = 0; k<V.numrows(); k++)
        {
            V[k][c]=vn[k][0]; //Subtracts orthogonalized col to V
            v[k][0]=V[k][c]; //v is redundant. more of a placeholder
        }
        if((c+1)<A.numcols())  //Matrix Out of Bounds Checker
        {
            for(int k = 0; k<A.numrows(); k++)
            {
                vn[k][0]=0;
                vc[k][0]=0;
                x[k][0]=A[k][c+1]; //Moves x onto next v
            }
        }
    }
    system("PAUSE");
    return V;
}

For testing purposes, I have been using the 2D Array: [[1,1,4],[1,4,2],[1,4,2],[1,1,0]]. Each column is its own 4x1 matrix. The matrices should be outputted as: [1,1,1,1]T, [-1.5,1.5,1.5,-1.5]T, and [2,0,0,-2]T respectively. What's happening now is that the first column comes out correctly (it's the same matrix), but the second and third come out to something that is potentially similar but not equal to their intended values.

Again, each time I call on the orthogonal method, it outputs something different. I think it's due to the numbers inputted in the proj() method, but I am not fully sure.

The apmatrix is from the AP college board, back when they taught cpp. It is similar to vectors or ArrayLists in Java.

Here is a link to apmatrix.cpp and to the documentation or conditions (probably more useful), apmatrix.h. Here is a link to the full code (I added visual markers to see what the computer is doing).

It's fair to assume that all custom methods work as intended (except maybe Matrix Regressions, but that's irrelevant). And be sure to enter the matrix using the enter method before trying to factorize. The code might be inefficient partly because I self-taught myself cpp not too long ago and I've been trying different ways to fix my code. Thank you for the help!


Solution

  • As said in comments:

    @AhmedFasih After doing more tests today, I have found that it is in-fact some >memory issue. I found that for some reason, if a variable or an apmatrix object >is declared within a loop, initialized, then that loop is reiterated, the >memory does not entirely wipe the value stored in that variable or object. This >is noted in two places in my code. For whatever reason, I had to set the >doubles a,b, and c to 0 in the proj method and apmatrixdh to 0 in the >mult method or they would store some value in the next iteration. Thank you so >much for you help!