Problem Statement

Any rectangular

m × n

{\displaystyle m\times n}

matrix

A

{\displaystyle {\mathit {A}}}

can be decomposed to a product of an orthogonal matrix

Q

{\displaystyle {\mathit {Q}}}

and an upper (right) triangular matrix

R

{\displaystyle {\mathit {R}}}

, as described in QR decomposition. Task Demonstrate the QR decomposition on the example matrix from the Wikipedia article: and the usage for linear least squares problems on the example from Polynomial regression. The method of Householder reflections should be used: Method Multiplying a given vector

a

{\displaystyle {\mathit {a}}}

, for example the first column of matrix

A

{\displaystyle {\mathit {A}}}

, with the Householder matrix

H

{\displaystyle {\mathit {H}}}

, which is given as reflects

a

{\displaystyle {\mathit {a}}}

about a plane given by its normal vector

u

{\displaystyle {\mathit {u}}}

. When the normal vector of the plane

u

{\displaystyle {\mathit {u}}}

is given as then the transformation reflects

a

{\displaystyle {\mathit {a}}}

onto the first standard basis vector which means that all entries but the first become zero. To avoid numerical cancellation errors, we should take the opposite sign of

a

1

{\displaystyle a_{1}}

: and normalize with respect to the first element: The equation for

H

{\displaystyle H}

thus becomes: or, in another form with Applying

H

{\displaystyle {\mathit {H}}}

on

a

{\displaystyle {\mathit {a}}}

then gives and applying

H

{\displaystyle {\mathit {H}}}

on the matrix

A

{\displaystyle {\mathit {A}}}

zeroes all subdiagonal elements of the first column: In the second step, the second column of

A

{\displaystyle {\mathit {A}}}

, we want to zero all elements but the first two, which means that we have to calculate

H

{\displaystyle {\mathit {H}}}

with the first column of the submatrix (denoted *), not on the whole second column of

A

{\displaystyle {\mathit {A}}}

. To get

H

2

{\displaystyle H_{2}}

, we then embed the new

H

{\displaystyle {\mathit {H}}}

into an

m × n

{\displaystyle m\times n}

identity: This is how we can, column by column, remove all subdiagonal elements of

A

{\displaystyle {\mathit {A}}}

and thus transform it into

R

{\displaystyle {\mathit {R}}}

. The product of all the Householder matrices

H

{\displaystyle {\mathit {H}}}

, for every column, in reverse order, will then yield the orthogonal matrix

Q

{\displaystyle {\mathit {Q}}}

. The QR decomposition should then be used to solve linear least squares (Multiple regression) problems

A

x

b

{\displaystyle {\mathit {A}}x=b}

by solving When

R

{\displaystyle {\mathit {R}}}

is not square, i.e.

m

n

{\displaystyle m>n}

we have to cut off the

m

− n

{\displaystyle {\mathit {m}}-n}

zero padded bottom rows. and the same for the RHS: Finally, solve the square upper triangular system by back substitution:

Let's start with the solution:

# R has QR decomposition built-in (using LAPACK or LINPACK)

a <- matrix(c(12, -51, 4, 6, 167, -68, -4, 24, -41), nrow=3, ncol=3, byrow=T)
d <- qr(a)
qr.Q(d)
qr.R(d)

# now fitting a polynomial
x <- 0:10
y <- 3*x^2 + 2*x + 1

# using QR decomposition directly
a <- cbind(1, x, x^2)
qr.coef(qr(a), y)

# using least squares
a <- cbind(x, x^2)
lsfit(a, y)$coefficients

# using a linear model
xx <- x*x
m <- lm(y ~ x + xx)
coef(m)