2017-02-20 3 views
2

J'écris un code pour la factorisation QR et, pour une raison quelconque, ma méthode orthogonale ne fonctionne pas comme prévu. Fondamentalement, ma méthode proj() produit des projections aléatoires. Voici le code:Orthogonalisation dans la factorisation QR produisant une matrice orthogonalisée légèrement inexacte

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; 
} 

J'ai testé la méthode par elle-même avec des entrées matricielles manuelles, et cela semble fonctionner correctement. Voici ma méthode orthogonale:

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; 
} 

Pour des fins de test, j'utilise le tableau 2D: [[1,1,4], [1,4,2], [1,4,2], [1,1,0]]. Chaque colonne est sa propre matrice 4x1. Les matrices doivent être sorties comme suit: [1,1,1,1] T, [-1,5,1,5,1,5, -1,5] T et [2,0,0, -2] T respectivement. Qu'est-ce qui se passe maintenant, c'est que la première colonne sort correctement (c'est la même matrice), mais la deuxième et la troisième sortent à quelque chose qui est potentiellement similaire, mais pas égal à leurs valeurs prévues.

Encore une fois, chaque fois que j'appelle la méthode orthogonale, elle produit quelque chose de différent. Je pense que c'est dû aux nombres entrés dans la méthode proj(), mais je ne suis pas entièrement sûr. L'apmatrice est issue du collège AP, à l'époque où elle enseignait le cpp. Il est similaire aux vecteurs ou aux ArrayLists en Java. Est un lien vers apmatrix.cpp et vers la documentation ou les conditions (probablement plus utiles), apmatrix.h. Here est un lien vers le code complet (j'ai ajouté des marqueurs visuels pour voir ce que fait l'ordinateur).

Il est juste de supposer que toutes les méthodes personnalisées fonctionnent comme prévu (à l'exception peut-être des régressions matricielles, mais ce n'est pas pertinent). Et assurez-vous d'entrer la matrice en utilisant la méthode enter avant d'essayer de factoriser. Le code peut être inefficace en partie parce que je me suis autodidacte cpp pas trop longtemps et j'ai essayé différentes façons de réparer mon code. Merci pour l'aide!

+0

Quelle est la différence entre les valeurs de ce qu'elles devraient être? –

+0

@AhmedFasih Chaque élément est toujours à + -1 de la valeur prévue. Une itération a retourné: [-1, 2, 2, -1] T pour la deuxième matrice et [2.8, -1, -1, 1.2] T. Parfois, la deuxième matrice est en fait les valeurs correctes. Notez également que la somme de tous les éléments de la deuxième matrice est égale à la somme de la matrice prévue. La troisième matrice serait aussi, sauf que deux des éléments sont supposés être 0. – Rickbox

+0

Ok, s'il ne s'agit pas d'erreurs d'arrondi numérique, alors je suspecterais un dépassement de tampon ou un problème de type d'accès mémoire puisque l'algorithme est entièrement déterministe mais montrant sorties aléatoires. Passez-le à travers valgrind ou certains désinfectants de Clang? –

Répondre

0

Comme dit dans les commentaires:

@AhmedFasih Après avoir fait plus d'essais aujourd'hui, je l'ai constaté qu'il est en fait une> question de la mémoire. J'ai trouvé que pour une raison quelconque, si une variable ou un objet apmatrix> est déclaré dans une boucle, initialisée, alors cette boucle est réitérée, la mémoire ne supprime pas entièrement la valeur stockée dans cette variable ou cet objet. Ceci est noté à deux endroits dans mon code. Pour une raison quelconque, j'ai dû mettre les doubles> a, b, et c à 0 dans la méthode proj et apmatrixdh à 0 dans la méthode> mult ou ils stockaient une valeur dans l'itération suivante. Merci beaucoup pour votre aide!