0

I am having trouble implementing the QR Decomposition using the Householder Reflection algorithm. I am getting the right numbers for the decomposition but the signs are incorrect. I can't seem to troubleshoot what's causing this issue since I can't find good pseudocode for the algorithm elsewhere on the internet. I have a feeling it has something to do with how I am choosing the sign for the tau variable, but I am following the pseudocode given to me from my textbook.

I have chosen to write this algorithm as a Fortran subroutine to prepare for an upcoming test on this material, and for now, it only computes the upper triangular matrix and not Q. Here is the output that this algorithm gives me:

  -18.7616634      -9.80723381      -15.7768545      -11.0864372    
   0.00000000      -9.99090481      -9.33576298      -7.53412533    
   0.00000000       0.00000000      -5.99453259      -9.80128098    
   0.00000000       0.00000000       0.00000000     -0.512612820

and here is the correct output:

   18.7616634       9.80723286       15.7768517       11.0864372    
   0.00000000       9.99090672       9.33576488       7.53412533    
   0.00000000       0.00000000       5.99453354       9.80128098    
   0.00000000       0.00000000       0.00000000      0.512614250

So it is basically just a sign issue it seems because the values themselves are otherwise correct.

! QR Decomposition with Householder reflections
! Takes a matrix square N x N matrix A
! upper triangular matrix R is stored in A
SUBROUTINE qrdcmp(A, N)
REAL, DIMENSION(10,10)  :: A, Q, I_N, uu
REAL gamma, tau, E
INTEGER N, & ! value of n for n x n matrix
        I,J,K ! loop indices

E = EPSILON(E) ! smallest value recognized by compiler for real

! Identity matrix initialization
DO I=1,N
    DO J=1,N
        IF(I == J) THEN
            I_N(I,J) = 1.0
        ELSE
            I_N(I,J) = 0.0
        END IF
    END DO
END DO


! Main loop
DO K=1,N

    IF(ABS(MAXVAL(A(K:N,K))) < E .and. ABS(MINVAL(A(K:N,K))) < E) THEN
        gamma = 0
    ELSE
        tau = 0
        DO J=K,N
            tau = tau + A(J,K)**2.0
        END DO
        tau = tau**(0.5) ! tau currently holds norm of the vector
        
        IF(A(K,K) < 0) tau = -tau
        
        A(K,K) = A(K,K) + tau
        gamma = A(K,K)/tau
        
        DO J=K+1,N
            A(J,K) = A(J,K)/A(K,K)
        END DO
        
        A(K,K) = 1
    END IF
    
    DO I=K,N
        DO J=K,N
            uu(I,J) = gamma*A(I,K)*A(J,K)
    END DO; END DO
    
    Q(K:N,K:N) = -1*I_N(K:N,K:N) + uu(K:N,K:N)
    
    A(K:N,K+1:N) = MATMUL(TRANSPOSE(Q(K:N,K:N)), A(K:N,K+1:N))
    A(K,K) = -tau
    
    
END DO

END SUBROUTINE qrdcmp

Any insight would be greatly appreciated. Also, I apologize if my code is not as easily readable as it could be.

Chris
  • 101

1 Answers1

2

In the construction of a reflector, you have the choice of two bisectors between the current column vector and the associated basis vector. Here this is done in the line IF(A(K,K) < 0) tau = -tau, as you already identified. "Unfortunately", the stable choice that avoids cancellation issues, usually also produces a sign flip in R as you observed. For instance, the decomposition of a matrix close to I will produce a result close to Q,R = -I,-I. If you want the R factor with positive diagonal, you have to shift the signs from R to Q in an extra step after completing the decomposition.

LAPACK had a short period where they used the choice producing the positive diagonal by default, this produced some strange errors/numerical stability problems where there previously where none.

Lutz Lehmann
  • 25,219
  • 2
  • 22
  • 51