IB Numerical Analysis - Linear least squares

9Linear least squares

IB Numerical Analysis

9 Linear least squares

Finally, we consider a slightly different kind of problem, whose solution involves

a different kind of factorization.

The question we are interested in is how we can “solve”

Ax = b,

where A ∈ R

m×n

, with m > n, and b ∈ R

and x ∈ R

Of course, this is an over-determined system. In general, this has no solution.

So we want to find the “best” solution to this system of equations. More precisely,

we want to find an

∗

∈ R

that minimizes

kAx − bk

. Note that we use the

Euclidean norm. This is the least squares problem.

Why would one be interested in this? Often in, say, statistics, we have some

linear model that says the outcome

is related to the input variables

, ··· , A

B = A

+ A

+ ···+ A

+ ε,

where the

are some unknown parameters, and

is some random fluctuation.

The idea is to do lots of experiments, say

of them, collect the data in

A ∈ R

m×n

and

b ∈ R

, and then find the

’s that predict these results the most accurately,

i.e. find the x such that Ax is as close to b as possible.

So how can we find this

∗

? First of all, we find a necessary and sufficient

condition for x

∗

to be a solution to this.

Theorem. A vector x

∗

∈ R

minimizes kAx − bk

if and only if

(Ax

∗

− b) = 0.

Proof. A solution, by definition, minimizes

f(x) = hAx − b, Ax − bi

= x

AAx − 2x

b + b

Then as a function of x, the partial derivatives of this must vanish. We have

∇f(x) = 2A

(Ax − b).

So a necessary condition is

(Ax − b).

Now suppose our

∗

satisfies

(

∗

− b

) = 0. Then for all

x ∈ R

, we write

x = x

∗

+ y, and then we have

kAx − bk

= kA(x

∗

+ y) −bk

= kAx

∗

− bk + 2y

(Ax − b) + kAyk

= kAx

∗

− bk + kAyk

≥ kAx

∗

− bk.

So x

∗

must minimize the Euclidean norm.

Corollary.

A ∈ R

m×n

is a full-rank matrix, then there is a unique solution

to the least squares problem.

Proof. We know all minimizers are solutions to

A)x = A

The matrix

being full rank means

y 6

0 ∈ R

implies

Ay 6

0 ∈ R

. Hence

A ∈ R

n×n

is positive definite (and in particular non-singular), since

Ax = (Ax)

(Ax) = kAxk

> 0

for x 6= 0. So we can invert A

A and find a unique solution x.

Now to find the minimizing x

∗

, we need to solve the normal equations

Ax = A

has full rank, then

is non-singular, and there is a unique solution. If

not, then the general theory of linear equations tells us there are either infinitely

many solution or no solutions. But for this particular form of equations, it turns

out we can show this is always consistent. So there will always be a solution.

However, solving these equations is not the recommended algorithm for finding

∗

for accuracy reasons. Instead the solution of these is based on orthogonal

matrices.

Definition

(Orthogonal matrix)

A square matrix

Q ∈ R

n×n

is orthogonal if

Q = I.

A key property of the orthogonal matrices is that we have

kQxk = kxk

for all

x ∈ R

. This is helpful since what we are trying to do in the least-squares

problem involves minimizing the Euclidean norm of something.

The idea is that for any orthogonal matrix, minimizing

kAx −bk

is the same

as minimizing

kQAx − Qbk

. So

should we pick? Recall that a simple kind of

matrices is triangular matrices. So we might want to get something, say, upper

triangular.

Definition

(QR factorization)

A QR factorization of an

m × n

matrix

is a

factorization of the form

A = QR,

where

Q ∈ R

m×m

is an orthogonal matrix, and

R ∈ R

m×n

is “upper triangular”

matrix, where “upper triangular” means

R =







··· R

0 R

··· R

0 0 ··· R

0 0 ··· 0







Note that every

A ∈ R

m×n

has a QR factorization, as we will soon show,

but this is not unique (e.g. we can multiply both Q and R by −1).

Once we have the QR factorization, we can multiply

kAx −bk

−1

and get an equivalent problem of minimizing

kRx −Q

. We will not go into

details, but it should be clear that this is not too hard to solve.

There is an alternative “skinny” QR factorization, with

, where we

remove redundant rows and columns such that

Q ∈ R

m×n

and

R ∈ R

n×n

. This

is good enough to solve the least squares problem.

We will see that if

has full rank, then the skinny QR is unique up to a

sign, i.e. if we require

> 0 for k = 1, ··· , n.

We shall demonstrate three standard algorithms for QR factorization:

(i) Gram-Schmidt factorization

(ii) Givens rotations

(iii) Householder reflections

The latter two are based on some geometric ideas, while the first is based on the

familiar Gram-Schmidt process.

Gram-Schmidt factorization

This targets the skinny version. So we stop writing the tildes above, and just

write Q ∈ R

m×n

and R ∈ R

n×n

As we all know, the Gram-Schmidt process orthogonalizes vectors. What we

are going to orthogonalize are the columns of A. We write

A =



··· a



, Q =



··· q



By definition of the QR factorization, we need

i=1

This is just done in the usual Gram-Schmidt way.

(i) To construct column 1, if a

6= 0, then we set

, R

= ka

Note that the only non-unique possibility here is the sign — we can let

−ka

and

−

instead. But if we require

0, then

this is fixed.

In the degenerate case,

, then we can just set

= 0, and the pick

any q

∈ R

with kq

k = 1.

(ii) For columns 1 < k ≤ n, for i = 1, ··· , k − 1, we set

= hq

, a

and set

= a

−

k−1

i=1

, a

If d

6= 0, then we set

and

= kd

In the case where

, we again set

= 0, and pick

to be anything

orthonormal to q

, ··· , q

k−1

In practice, a slightly different algorithm (modified Gram-Schmidt process)

is used, which is (much) superior with inexact arithmetic. The modified Gram-

Schmidt process is in fact the same algorithm, but performed in a different order

in order to minimize errors.

However, this is often not an ideal algorithm for large matrices, since there

are many divisions and normalizations involved in computing the

, and the

accumulation of errors will cause the resulting matrix Q to lose orthogonality.

Givens rotations

This works with the full QR factorization.

Recall that in R

, a clockwise rotation of θ ∈ [−π, π] is performed by



cos θ sin θ

−sin θ cos θ







α cos θ + β sin θ

−α sin θ + β cos θ



By choosing θ such that

cos θ =

+ β

, sin θ =

+ β

this then becomes



+ β



Of course, by choosing a slightly different

, we can make the result zero in the

first component and

+ β

in the second.

Definition

(Givens rotation)

, where

m >

2, we define the Givens

rotation on 3 parameters 1 ≤ p < q ≤ m, θ ∈ [−π, π] by

Ω

[p,q]







cos θ sin θ

−sin θ cos θ







∈ R

m×m

where the sin and cos appear at the p, qth rows and columns.

Note that for

y ∈ R

, the effect of Ω

[p,q]

only alters the

and

components.

In general, for B ∈ R

m×n

, then Ω

[p,q]

B only alters the p and q rows of B.

Moreover, just like the

case, given a particular

z ∈ R

, we can choose

such that the qth component (Ω

[p,q]

= 0.

Hence,

A ∈ R

m×n

can be transformed into an “upper triangular” form by

applying

mn −

n(n+1)

Givens rotations, since we need to introduce

many

zeros. Then

···Q

A = R.

How exactly do we do this? Instead of writing down a general algorithm, we

illustrate it with a matrix A ∈ R

4×3

. The resultant R would look like

R =







0 R

0 0 R

0 0 0







We will apply the Givens rotations in the following order:

Ω

[3,4]

Ω

[3,4]

Ω

[2,3]

Ω

[1,4]

Ω

[1,3]

Ω

[1,2]

A = R.

The matrix A transforms as follows:

× × ×

Ω

[1,2]

× × ×

× ×

× × ×

Ω

[1,3]

× × ×

× ×

× × ×

Ω

[1,3]

× × ×

× ×

Ω

[2,3]

× × ×

× ×

0 0

× ×

Ω

[3,4]

× × ×

× ×

0 0

Ω

[3,4]

× × ×

× ×

0 0

0 0 0

Note that when applying, say, Ω

[2,3]

, the zeros of the first column get preserved,

since Ω

[2,3]

only mixes together things in row 2 and 3, both of which are zero in

the first column. So we are safe.

Note that this gives us something of the form

···Q

A = R.

We can obtain a QR factorization by inverting to get

A = Q

···Q

However, we don’t really need to do this if we just want to do use this to solve

the least squares problem, since to do so, we need to multiply by

, not

which is exactly Q

···Q

Householder reflections

Definition

(Householder reflection)

For

u 6

0 ∈ R

, we define the House-

holder reflection by

= I − 2

∈ R

m×m

Note that this is symmetric, and we can see

. So this is indeed

orthogonal.

To show this is a reflection, suppose we resolve

as the perpendicular and

parallel parts as x = αu + w ∈ R

, where

α =

, u

w = 0.

Then we have

x = −αu + w.

So this is a reflection in the m −1-dimensional hyperplane u

y = 0.

What is the cost of computing H

z? This is evaluated as

z = z − 2

This only requires O(m) operations, which is nice.

Proposition.

A matrix

A ∈ R

m×n

can be transformed into upper-triangular

form by applying n Householder reflections, namely

···H

A = R,

where each

introduces zero into column

and leaves the other zeroes alone.

This is much less multiplications needed than the number needed for Givens

rotations, in general.

To show this, we first prove some lemmas.

Lemma.

Let

a, b ∈ R

, with

a 6

, but

kak

kbk

. Then if we pick

a −b

then

a = b.

This is obvious if we draw some pictures in low dimensions.

Proof. We just do it:

a = a −

2(kak − a

kak

− 2a

b + kbk

(a − b) = a − (a − b) = b,

where we used the fact that kak = kbk.

We will keep applying the lemma. Since we want to get an upper triangular

matrix, it makes sense to pick

to have many zeroes, and the best option would

be to pick b to be a unit vector.

So we begin our algorithm: we want to clear the first column of

. We let

be the first column of

, and assume

a ∈ R

is not already in the correct form,

i.e. a is not a multiple of e

. Then we define

u = a ∓ kake

where either choice of the sign is pure-mathematically valid. However, we will

later see that there is one choice that is better when we have to deal with inexact

arithmetic. Then we end up with

a = H

a = ±kake

Now we have

A =

× × ×

× ×

To do the next step, we need to be careful, since we don’t want to destroy the

previously created zeroes.

Lemma. If the first k − 1 components of u are zero, then

(i) For every x ∈ R

, H

x does not alter the first k −1 components of x.

(ii) If the last (m − k + 1) components of y ∈ R

are zero, then H

y = y.

These are all obvious from definition. All these say is that reflections

don’t affect components perpendicular to

, and in particular fixes all vectors

perpendicular to u.

Lemma. Let a, b ∈ R

, with

























but

j=k

Suppose we pick

u = (0, 0, ··· , 0, a

− b

, ··· , a

− b

)

Then we have

a = (a

, ··· , a

k−1

, ··· , b

This is a generalization of what we’ve had before for

= 1. Again, the proof

is just straightforward verification.

Now we can mess with the second column of

. We let

be the second

column of

, and assume

, ··· , a

are not all zero, i.e. (0

, a

, ··· , a

)

not a multiple of e

. We choose

u = (0, a

∓ γ, a

, ··· , a

)

where

γ =

j=2

Then

a = (a

, ±γ, 0, ···).

Also, by the previous lemma, this does not affect anything in the first column

and the first row. Now we have

A =

× × ×

0 0

Suppose we have reached

k−1

···H

, where the first

k −

1 rows are of the

correct form, i.e.

k−1

···H

A =

k − 1

We consider the

th column of

k−1

···H

, say

. We assume

(0, ··· , 0, a

, ··· , a

)

is not a multiple of e

. Choosing

u = γ(0, ··· , 0, a

∓ γ, a

k+1

, ··· , a

)

, γ =

j=k

we find that

a = (a

, ··· , a

k−1

, ±γ, 0, ··· , 0)

Now we have

···H

A =

and H

did not alter the first k −1 rows and columns of H

k−1

···H

There is one thing we have to decide on — which sign to pick. As mentioned,

these do not matter in pure mathematics, but with inexact arithmetic, we

should pick the sign in

∓ γ

such that

∓ γ

has maximum magnitude, i.e.

∓ γ

sgn

(

)

. It takes some analysis to justify why this is the right

choice, but it is not too surprising that there is some choice that is better.

So how do Householder and Givens compare? We notice that the Givens

method generates zeros at one entry at a time, while Householder does it column

by column. So in general, the Householder method is superior.

However, for certain matrices with special structures, we might need the

extra delicacy of introducing a zero one at a time. For example, if

is a band

matrix, then it might be beneficial to just remove the few non-zero entries one

at a time.