Finding vector $x$ so that $Ax=b$ using Householder reflections.

Question

Assume $n\times m$ matrix $A$ and vector $b$ are given. I am looking for $x$ that satisfies $Ax=b$ in terms of linear least squares problem. Let $A=\begin{bmatrix} 1 & 1 & 1 \\0 & 1 & 1\\0 & 0 & 1\end{bmatrix}$ and $b=\begin{bmatrix} 3\\2\\1\end{bmatrix}$. Using Householder reflections on both matrices columns and then calcualting $x$, well, doesn't work. I know $A$ is uppertriangular here, but does that mean there are there some conditions under which I can't use Householder method? I wrote a code in Octave that deals with above-mentioned task, but it won't always work properly :/ That's why I'm asking.

Or is it something wrong with the code itself?

function [x,R,B]=householder(A,b)
n=columns(A);
B=P=R=x=[];
d=b;
for i=[1:1:n]
    u=A(:,i);
    v=u;            %Householder method in use
    if (i>1)        %starting form 2nd column, I make set all values above i-th to 0
        v([1:i-1],:)=0;
    endif
    v([i],:)=v([i])+norm(v);
    B=[B,v];                %don't mind B
    P=[P,(u-v)];
    d=d-2/(norm(v))^2*v'*d*v    %I transform the result vector
endfor  
R=P([1:n],[1:n]);
d=d([1:n]);
for i=[n:-1:1]
    y=d([i])/R([i],[i]);
    if (i>1)
        for j=[i-1:-1:1]        %I calculate x by manipulating d
            d([j])=d([j])-R([j],[i])*y;
        endfor
    endif
    x=[y;x];
endfor

Heeeelp.

Answer 1

Try this instead (for MATLAB but it should work in Octave too):

function [x, R, c] = hh(A, b)

[m, n] = size(A);
assert(m >= n);
assert(size(b, 1) == m);

R = A;
c = b;

for i = 1 : n
    v = R(i:end,i);
    alpha = -norm(v);
    if (v(1) < 0), alpha = -alpha; end
    v(1) = v(1) - alpha;
    v = v / norm(v);
    R(i,i) = alpha; R(i+1:end,i) = 0;
    R(i:end,i+1:end) = R(i:end,i+1:end) ...
                       - 2 * v * (v' * R(i:end,i+1:end));
    c(i:end,:)         = c(i:end,:) ...
                       - 2 * v * (v' * c(i:end,:));
end

dR = abs(diag(R));
if (min(dR) / max(dR) < eps * max(m, n))
    warning('Matrix A is probably badly scaled or ill-conditioned.');
end

x = c(1:n,:);
for i = n : -1 : 1
    for j = i + 1 : n
        x(i,:) = x(i,:) - R(i,j) * x(j,:);
    end
    x(i,:) = x(i,:) / R(i,i);
end

The function assumes $A$ to be $m\times n$ with $m\geq n$ and $b$ to be $m\times p$ and computes the solution of $Ax=b$ (in the least squares sense if $Ax=b$ is not consistent) provided that the rank of $A$ is $n$.