Lecture 10: Stabilized Orthogonal Factorizations

1.Conditioning of linear algebra problems

Recall that the relative condition number of a mathematical problem $f:X\to Y$ characterizes the amplification by $f$ of perturbations in its argument

Linear combination. The basic operation of linear combination $𝑨𝒙$, $𝑨\in {ℂ}^{m\times n}$, seen as the problem ${ℂ}^n\stackrel{𝒇}{\to }{ℂ}^m$ has the condition number

The matrix norm property $||𝑨𝒚||⩽||𝑨||||𝒚||$ can be used to obtain

leading to

where $\kappa \left(𝑨\right)$ is the condition number of the matrix $𝑨$. If $𝑨$ is of full rank with $m>n$, the 2-norm condition number is given by the ratio of largest to smallest singular values.

By convention, if $𝑨$ is singular, the condition number $\kappa \left(𝑨\right)=\infty$.

Coordinate transformation. For $𝑨\in {ℂ}^{m\times m}$ of full rank, the coordinates of vector $𝒃\in {ℂ}^m$ expressed in the $𝑰$ basis can be transformed its coordinates $𝒙\in {ℂ}^m$ in the $𝑨$ basis by solving the linear system $𝑨𝒙=𝑰𝒃$, with the solution $𝒙={𝑨}^{-1}𝒃$ (so written formally, even though the inverse is almost never explicitly computed). This is simply another linear combination of the columns of ${𝑨}^{-1}$, hence the problem $𝒇:{ℂ}^m\to C^m$, $𝒇\left(𝒃\right)={𝑨}^{-1}𝒃$ again has a condition number bounded by the condition number of the matrix $𝑨$.

Operator perturbation. Instead of changing the input data as above, the linear mapping represented by the matrix $𝑨\in {ℂ}^{m\times n}$ might itself be perturbed. Two mathematical problems may now be formulated:

1. For fixed $𝒃\in {ℂ}^m,𝒇:{ℂ}^{m\times n}\to {ℂ}^n,𝒇(𝑨,𝒃)={𝑨}^{+}𝒃=𝒙$. Perturbation of the input $𝑨$ induces perturbation of $𝒙$ in order for $𝒃$ to be kept fixed

   ${\displaystyle (𝑨+\delta 𝑨)(𝒙+\delta 𝒙)=𝒃.}$

   Using $𝑨𝒙=𝒃$, and assuming that $\delta 𝑨\delta 𝒙$ is negligible gives

   ${\displaystyle 𝑨\delta 𝒙+\delta 𝑨𝒙=𝟎\Rightarrow \delta 𝒙=-{𝑨}^{+}\delta 𝑨𝒙,}$

   hence the relative condition number is

   ${\displaystyle \kappa =\frac{||{𝑨}^{+}\delta 𝑨𝒙||}{||𝒙||}\cdot \frac{||𝑨||}{||\delta 𝑨||}⩽\frac{||{𝑨}^{+}||||\delta 𝑨𝒙||}{||𝒙||}\cdot \frac{||𝑨||}{||\delta 𝑨||}⩽\frac{||{𝑨}^{+}||||\delta 𝑨||||𝒙||}{||𝒙||}\cdot \frac{||𝑨||}{||\delta 𝑨||}=\kappa \left(𝑨\right).}$

For all above linear algebra problems the condition number is bounded by the associated matrix condition number. Unitary matrices $𝑸\in {ℂ}^{m\times m}$ have $\kappa \left(𝑸\right)=1$, and define an orthonormal basis for ${ℂ}^m$. A rank-deficient matrix $𝒁\in {ℂ}^{m\times m}$ has $\kappa \left(𝒁\right)=\infty$, and corresponds to a linearly dependent vector set $\{{𝒛}_1,\dots ,{𝒛}_m\}$. The behavior of many numerical approximation procedures based upon linear combinations is determined by condition number of the basis set.

∴ 

function Vandermonde(a,b,m)
  t=LinRange(a,b,m); v=ones(m,1); V=copy(v)
  for j=2:m
    v = v .* t; V=[V v]
  end
  return V
end;

∴ 

∴ 

function VandermondeC(m)
  δ=π/(2*m); ϴ=LinRange(δ,π-δ,m)
  t=cos.(ϴ)
  v=ones(m,1); V=copy(v)
  for j=2:m
    v = v .* t; V=[V v]
  end
  return V
end;

∴ 

∴ 

function Triangular(a,b,m)
  x=LinRange(a,b,m); T=ones(m,1); Tj=copy(T); t=copy(x)
  for j=2:m
    Tj = Tj .* (t .- x[j-1]); T=[T Tj]
  end
  return T
end;

∴ 

∴ 

function TriangularC(m)
  δ=π/(2*m); ϴ=LinRange(δ,π-δ,m)
  x=cos.(ϴ); T=ones(m,1); Tj=copy(T); t=copy(x)
  for j=2:m
    Tj = Tj .* (t .- x[j-1]); T=[T Tj]
  end
  return T
end;

∴ 

$•$

Figure 1. Monomial basis with: (o) uniform sampling, (x) Chebyshev sampling. Triangular basis with: (+) uniform sampling, (*) Chebyshev sampling.

∴	mr=5:5:100; κVDMU=log10.(cond.(Vandermonde.(-1,1,mr)));

∴	κVDMC=log10.(cond.(VandermondeC.(mr)));

∴	κTU=log10.(cond.(Triangular.(-1,1,mr)));

∴	κTC=log10.(cond.(TriangularC.(mr)));

∴

∴	x=collect(mr); clf();

∴	plot(x,κVDMU,"o-",x,κVDMC,"x-",κTU,"+-",κTC,"*-");

∴	grid("on"); title("Condition number κ of polynomial bases");

∴	xlabel("Number of sample points"); ylabel("lg(κ)");

∴	pre=homedir()*"/courses/MATH661/images/";

∴	savefig(pre*"PolyBasesCondNr.eps");

∴

2.Orthogonal factorization through Householder reflectors

The Gram-Schmidt procedure constructs an orthogonal factorization by linear combinations of the column vectors of $𝑨\in {ℂ}^{m\times n}$, $m⩾n$, $\mathrm{rank}\left(𝑨\right)=n$

In exact arithmetic $C\left(𝑸\right)=C\left(𝑨\right)$ by construction, and $\kappa \left(𝑸\right)=1$, but the sequence of multiplications with ${𝑹}_1,\dots ,{𝑹}_n$ might amplify perturbations in $𝑨$ (due for example to floating point representation errors or inherent uncertainty in knowledge of $𝑨$). The problem $𝒇:{ℂ}^{m\times n}\to C^{m\times n}\times {ℂ}^{n\times n},$ $𝑨\stackrel{𝒇}{\to }𝑸,𝑹$ has condition number

and numerical experimentation (Fig. 2) readily exhibits large condition numbers.

An alternative approach is to obtain an orthogonal factorization through unitary transformations

Unitary transformations do not modify vector 2-norms or relative orientations

and are hence said to be isometric. In Euclidean space reflections and rotations are isometric.

$\circ$

Figure 2. $QR$-conditioning: (o) modified Gram-Schmidt, (x) Householder.

Estimation of the $QR$-condition number by numerical experimentation: generate $N$ random perturbations of a matrix, compute the factorization of each perturbed matrix, and choose the maximum encountered value as $\kappa$.

∴	function mgs(A) m,n=size(A); Q=copy(A); R=zeros(n,n) for i=1:n R[i,i]=sqrt(Q[:,i]'*Q[:,i]) if (R[i,i]<eps()) break end Q[:,i]=Q[:,i]/R[i,i] for j=i+1:n R[i,j]=Q[:,i]'*A[:,j] Q[:,j]=Q[:,j]-R[i,j]*Q[:,i] end end return Q,R end;

∴	function QRcondGS(N,A,ε) Q,R=mgs(A); κ=1; normA=norm(A); normR=norm(R) for k=1:N δA=ε*randn(size(A)) δF=qr(A+δA); δQ=Array(δF.Q)-Q; δR=Array(δF.R)-R κ = max(κ, (norm(δQ)+norm(δR)/normR)*normA/norm(δA) ) end return κ end;

∴	mr=10:10:100; κQRGS = log10.(QRcondGS.(100,randn.(mr,mr),1.0e-6));

∴	x=collect(mr); clf(); plot(x,κQRGS,"o");

∴

∴	function QRcond(N,A,ε) F=qr(A); Q=Array(F.Q); R=Array(F.R) κ=1; normA=norm(A); normR=norm(R) for k=1:N δA=ε*randn(size(A)) δF=qr(A+δA); δQ=Array(δF.Q)-Q; δR=Array(δF.R)-R κ = max(κ, (norm(δQ)+norm(δR)/normR)*normA/norm(δA) ) end return κ end;

∴	κQR = log10.(QRcond.(100,randn.(mr,mr),1.0e-6));

∴	x=collect(mr); plot(x,κQR,"x");

∴	grid("on"); title("Condition number of QR-factorization of A");

∴	xlabel("Dimension of A"); ylabel("log10(κ)");

∴	savefig(pre*"QRcond.eps");

∴

Construction of an isometric reflection transformation suitable for a $QR$ factorization is represented in Fig. 3. Let vector $𝒙\in {ℂ}^{m+1-k}$ represent the portion of the $k^{\mathrm{th}}$ column from the diagonal downwards in stage $k$ of reduction of $𝑨\in {ℂ}^{m\times n}$ to upper triangular form

The next stage of in reduction to upper triangular form is the isometric transformation of $𝒙$ into $\pm ||𝒙||{𝒆}_1$, with ${𝒆}_1\in {ℂ}^{m+1-k}$ the unit vector along the first direction. With $𝒗=\pm ||𝒙{𝒆}_1||-𝒙$, $𝒒=𝒗/||𝒗||$, the projection of $𝒙$ onto the span of $𝒗$, $C\left(𝒗\right)$ is

and its complementary projector onto $N\left({𝒗}^{\ast }\right)$ is

The reflector transforming $𝒙$ into $\pm ||𝒙||{𝒆}_1$ is obtained by doubling the above displacements, and is known as a Householder reflector

Of the two possibilities $\pm ||𝒙||{𝒆}_1$, the choice

avoids loss of floating accuracy $𝒙\cong ||𝒙||{𝒆}_1$. For $𝒙\in {ℂ}^{m+1-k}$, $\mathrm{sign}\left(x_1\right)=\exp \left(\arg \left(x_1\right)\right)$.

Figure 3. Geometry of Householder reflector

The resulting Householder $QR$-factorization is given

∴ 

function HouseholderQR(A)
  m,n=size(A)
  Q=zeros(m,n); R=copy(A)
  for k=1:n
    x=R[k:m,k]
    e1=zeros(size(x)); e1[1]=1
    v=sign(x[1])*norm(x)*e1+x
    q=v/norm(v); Q[k:m,k]=q
    for j=k:n
      aj=R[k:m,j]; c=2*q'*aj
      R[k:m,j]=aj.-c*q
    end
  end
  return Q,R
end;

∴ 

Note that the above implementation does not return the $𝑸$ matrix, but rather the ${𝑸}_1,\dots ,{𝑸}_n$ reflectors from which $𝑸$ can be reconstructed if needed. Usually though, the $𝑸$ matrix itself is not required, but rather the product $𝑸𝒖$ which can readily be evaluated as ${𝑸}_n\dots {𝑸}_1𝒖$. The Householder reflector algorithm is typically the default procedure in $QR$-factorizations implemented in software systems, and as seen in (Fig. 2), leads to much better conditioning.

3.Orthogonal factorization through Given rotators

An alternative approach to orthogonal factorization utilizes isometric rotation transformations of the form

with the rotation angle $\theta$ chosen to nullify the subdiagonal element $(i,k)$, $i>k$

Composition of repeated rotations ${𝑸}_{ik}=𝑹(i,k,{\theta }_{ik})$ can be organized to lead to an upper triangular matrix

Whereas Householder reflectors work on entire vectors, Givens rotators nullify individual subdiagonal elements. For full matrices Householder reflectors typically require fewer floating point operations, but the special structure of a sparse matrix is better exploited by use of Givens rotators.

∴ 

function GivensQR(A)
  m,n=size(A)
  Q=zeros(m,n); R=copy(A)
  for k=1:n
    for i=k+1:m
      θ = atan(-R[i,k],R[k,k]); Q[i,k]=
      c = cos(θ); s = sin(θ)
      for j=k:n
        u = R[k,j]; v = R[i,j]
        R[k,j]=c*u-s*v
        R[i,j]=s*u+c*v
      end
    end
  end
  return Q,R
end;

∴ 

As in the Householder implementation the above implementation returns data to reconstruct $𝑸$ if needed.