Problem implementing a QR factorization

Question

I'm trying to write a Fortran subroutine to compute a QR factorization using the Householder method. To test my routine, I compute the factorization of the following matrix: $$ A = \begin{pmatrix} 12 & -51 & 4 \\ 6 & 167 & -68 \\ -4 & 24 & -41 \end{pmatrix}, $$ which, if done correctly, will reduce to the following upper triangular matrix: $$ R = \begin{pmatrix} 14 & 21 & -14 \\ 0 & 175 & -70 \\ 0 & 0 & 35 \end{pmatrix}. $$ However, the matrix I actually get is: $$ R = \begin{pmatrix} -14 & -21 & 14 \\ 0 & -175 & 70 \\ -0 & 0 & 35 \end{pmatrix}, $$ which looks almost correct, except for some strange sign changes. I've been staring at my subroutine all day trying to see where these sign changes are being introduced, but I can't identify the problem.

My algorithm is as follows:

$ for \:\: k \:=\: 1\: to\: n $

$ \qquad x(k:m) = A(k:m,k) $

$ \qquad v(k:m) = \mathtt{sign}(x(k))||x(k:m)||_{2}e1 + x(k:m) $

$ \qquad v(k:m) = v(k:m)/||v(k:m)||_{2} $

$ \qquad A(k:m,k:n) = A(k:m,k:n) - 2vv^{\top}A(k:m,k:n) $

To calculate the factor $2vv^{\top}A(k:m,k:n),$ I made another subroutine called outer_product to compute the outer product of $v$ with itself, i.e. $vv^{\top}$, and then matrix multiply the result into my submatrix $A(k:m,k:n)$. However, I'm not sure if this is legitimate - I suspect herein lies the problem.

I would really appreciate it if someone could glance at my code to see if there is any obvious reason for the incorrect sign changes:

integer, parameter :: dp = selected_real_kind(15)

integer, intent(in) :: m, n
real(dp), dimension(m,n), intent(inout) :: A
real(dp), dimension(m,m), intent(out) :: Q

integer :: k
real(dp) :: two_norm
real(dp), dimension(m) :: x, e1
real(dp), dimension(m,n) :: v
real(dp), dimension(m,m) :: outprod_vv

v = 0.0_dp

do k=1,m
    Q(k,k) = 1.0_dp
end do

!Householder triangularization
do k=1,n

    e1(k) = 1.0_dp

    x(k:m) = A(k:m,k)
    v(k:m,k) = sign( sqrt(dot_product(x(k:m),x(k:m))), x(k) )* &
        e1(k:m) + x(k:m)

    v(k:m,k) = v(k:m,k)/(sqrt(dot_product(v(k:m,k),v(k:m,k))))
    call outer_product(v(k:m,k), m-k+1, outprod_vv(k:m,k:m))

    A(k:m,k:n) = A(k:m,k:n) - &
        2.0_dp*matmul(outprod_vv(k:m,k:m), A(k:m,k:n)) 

    !Form Q implicitly    
    Q(k:m,k:m) = Q(k:m,k:m) - 2.0_dp* &
        matmul(outprod_vv(k:m,k:m), Q(k:m,k:m))

end do

Q = transpose(Q)

Answer 1

A quick glance tells me your routine is fine. But more importantly, since QR factorization is not unique, the $R$ your routine outputs is perfectly fine-- it just means that the $Q$ has absorbed a negative. Also, for comparison, the Householder routine in Matlab gives the same $R$ as your Fortran routine.

If you implement the modified Gram-Schmidt (MGS) routine, you will obtain the $R$ you originally sought, since MGS sets the diagonal entries to positive values.

In practice, this difference in sign doesn't make a difference, since one usually cares more about the product $QR$ than about the individual factors.

Answer 2

A Householder reflection $Q:=I-2vv^T$ which transforms a vector $x$ to a distinct vector $y$ (such that $\|x\|_2=\|y\|_2$) can be constructed by setting $v:=(x-y)/\|x-y\|_2$. We have $$ Qx=\left[1-\frac{2(x-y)^Tx}{(x-y)^T(x-y)}\right]x+\frac{2(x-y)^Tx}{(x-y)^T(x-y)}y. $$ Using the above assumptions, it is straightforward to verify that $Qx=y$ (the coefficients by $x$ and by $y$ are, respectively, zero and one).

So, to zero out all the components of $x$, $y:=\sigma e_1$, where $|\sigma|=\|x\|_2$. Since one usually wants to avoid severe cancelation in computing $v$, $\sigma$ is chosen so that the sign of the first component of $x$ is the same as the sign of $\sigma$. Note that, in $Qx=\sigma e_1$, $\sigma$ is negative if the first component of $x$ is positive so that explains why HHQR creates the upper triangular matrix with some negative diagonal entries.

If you would like to have a positive diagonal in $R$, consider multiplying by $-Q=2vv^T-I$ instead of by $Q$.

Since you showed the code, I would have some remarks:

• You do not need a special vector for $e_1$ since adding its multiple to a vector does nothing but modifies one entry.

• To compute $QB=B-2vv^TB$, it is better to implement it as $QB=B-2v(v^TB)$ instead of $QB=B-2(vv^T)B$. I believe you can see why if, e.g., $m\gg n$.

• Although in HHQR, it is better to do not store $Q$ explicitly (think again about the case $m\gg n$; not mentioning numerical reasons now), your $Q$ is certainly incorrect. You should update all the columns of $Q$ instead of $k$ to $m$.

• Since by construction of $Q$, you know what should be the content of $A(k:m,k)$ (its first component is up to sign equal to the norm of $v$ and the rest is zero). So computing it explicitly by applying $Q$ is, firstly, waste of time, and, secondly, instead of having zeros in rows $k+1,\ldots,m$, you are very likely to get some small junk due to roundoff in the strictly lower triangular part of $R$.

For inspiration, you can have a look on my answer here; you might be intereseted there in particular in the implementation of the orthogonalization loop. I would also suggest to have a look on Chapter 5 here.