Orthogonality

Example

Remember the example

A = [a_{1}, a_{2}, a_{3}, a_{4}] = [ \begin{matrix}1 & 2 & 3 & 4 \\ 4 & 3 & 2 & 1 \\ −10 & −5 & 0 & 5 \end{matrix} ] , b = [ \begin{matrix}1 \\ 1 \\ 1 \end{matrix} ] , c = [ \begin{matrix}1 \\ 1 \\ −1 \end{matrix} ] .

The Gauss elimination method on [A|b c] (in Matlab: rref([A b c])) gives

[Â|b̂ ĉ] = [ \begin{matrix}1 & 0 & −1 & −2 \\ 0 & 1 & 2 & 3 \\ 0 & 0 & 0 & 0 \end{matrix} | \begin{matrix}0 & −1/5 \\ 0 & 3/5 \\ 1 & 0 \end{matrix} ] .

We conclude that {a_{1}, a_{2}} are linearly independent and that they are a basis for the linear space

V = R(A) = S(a_{1}, a_{2}, a_{3}, a_{4}) = S(a_{1}, a_{2}), dim(V) = 2.

We also see that b ∉ V because Ax = b has no solution, and that c ∈ V because Ax = c has the multiple solutions

x = \underbrace{s [ \begin{matrix}1 \\ −2 \\ 1 \\ 0 \end{matrix} ] + t [ \begin{matrix}2 \\ −3 \\ 0 \\ 1 \end{matrix} ]}_{=x_{h}∈N(A)} + \underbrace{[ \begin{matrix}−1/5 \\ 3/5 \\ 0 \\ 0 \end{matrix} ]}_{=x_{p}} = x_{h} + x_{p}.

We write the basis {a_{1}, a_{2}} as a matrix

B = [a_{1}, a_{2}] = [ \begin{matrix}1 & 2 \\ 4 & 3 \\ −10 & −5 \end{matrix} ] .

and note that V = R(A) = R(B) = S(a_{1}, a_{2}) is a two-dimensional subspace of ℝ^{3} (a plane through the origin).

Figure is missing!!

In order to express c as a linear combination of the basis, i.e.,

c = Bx = x_{1}a_{1} + x_{2}a_{2},

we solve the equation Bx = c and get the unique solution x = [ \begin{matrix}−1/5 \\ 3/5 \end{matrix} ], so that

c = −\frac{1}{5} [ \begin{matrix}1 \\ 4 \\ −10 \end{matrix} ] + \frac{3}{5} [ \begin{matrix}2 \\ 3 \\ −5 \end{matrix} ] = [ \begin{matrix}1 \\ 1 \\ −1 \end{matrix} ] .

Ortogonal projektion på V

(AMBS 42.35–38) I allmänhet betraktar vi en bas B = [a_{1}, …, a_{n}], m × n, för ett underrum V = R(B) = S(a_{1}, …, a_{n}) till ℝ^{m}, n ≤ m, och en vektor b ∈ ℝ^{m}, b ∉ V. I exemplet: m = 3, n = 2.

Vi vill hitta en vektor P b ∈ V sådan att b − P b ⊥ V, dvs

(b − P b, v) = 0 ∀v ∈ V.  (1)

Det räcker att testa ortogonalitet mot alla basvektorer, dvs

(b − P b, a_{j}) = 0, j = 1, …, n.  (2)

För om (2) gäller så skriver vi v = ∑_{j=1}^{n} α_{j}a_{j} och

(b − P b, ∑_{j=1}^{n} α_{j}a_{j}) = ∑_{j=1}^{n} α_{j} \underbrace{(b − P b, a_{j})}_{=0} = 0,

så att (1) gäller. Det är klart att (1) medför (2).

Nu beräknar vi P b. Vi gör en ansats: P b = ∑_{i=1}^{n} x_{i}a_{i} = Bx. Vi måste bestämma koefficienterna x_{i}. Ekvation (2) ger

\begin{matrix}(P b, a_{j}) = (b, a_{j}), j = 1, …, n \\ ∑_{i=1}^{n} x_{i}(a_{i}, a_{j}) = (b, a_{j}), j = 1, …, n. \end{matrix}

Detta är ett ekvationsystem för x_{i}. Alternativt kan man skriva detta på matrisform

(Bx, a_{j}) = (b, a_{j}), j = 1, …, n,
a_{j}^{T} Bx = a_{j}^{T} b, j = 1, …, n,
B^{T} Bx = B^{T} b.

Detta är ett system av n ekvationer i n obekanta med symmetrisk matris: B^{T} B.

\underbrace{\underbrace{B^{T}}_{n×m} \underbrace{B}_{m×n}}_{n×n} \underbrace{x}_{n×1} = \underbrace{\underbrace{B^{T}}_{n×m} \underbrace{b}_{m×1}}_{n×1}

I exemplet:

B^{T} B = [ \begin{matrix}1 & 4 & −10 \\ 2 & 3 & −5 \end{matrix} ] [ \begin{matrix}1 & 2 \\ 4 & 3 \\ −10 & −5 \end{matrix} ] = [ \begin{matrix}117 & 64 \\ 64 & 38 \end{matrix} ]

B^{T} b = [ \begin{matrix}1 & 4 & −10 \\ 2 & 3 & −5 \end{matrix} ] [ \begin{matrix}1 \\ 1 \\ 1 \end{matrix} ] = [ \begin{matrix}−5 \\ 0 \end{matrix} ]

Ekvationssystemet B^{T} Bx = B^{T} b blir

[ \begin{matrix}117 & 64 \\ 64 & 38 \end{matrix} ] [ \begin{matrix}x_{1} \\ x_{2} \end{matrix} ] = [ \begin{matrix}−5 \\ 0 \end{matrix} ]

och leder till trappstegsmatrisen [ \begin{matrix}117 & 64 \\ 0 & \frac{350}{117} \end{matrix} | −5_{\frac{320}{117}} ]. Vi får unik lösning

\{ \begin{matrix}x̂_{2} = \frac{320}{350} ≈ 0.9143 \\ x̂_{1} ≈ −0.5429 \end{matrix}

dvs x̂ = [ \begin{matrix}−0.543 \\ 0.914 \end{matrix} ] och P b = Bx̂ ≈ [ \begin{matrix}1.286 \\ 0.571 \\ 0.857 \end{matrix} ]

Vi går tillbaka till det allmänna fallet. Vi vill lösa B^{T} Bx = B^{T} b. Vi ska visa att matrisen B^{T} B är icke-singulär. Eftersom den är en kvadratisk matris, n × n, räcker det att visa att N(B^{T} B) = {0} enligt Sats 42.6.

Vi kollar detta genom att lösa det homogena systemet B^{T} Bx = 0. Multiplicera med x^{T}:

x^{T} B^{T} Bx = 0
(Bx)^{T} (Byx) = 0
∥Bx∥^{2} = 0
Bx = 0
x = 0 ty a_{1}, …, a_{n} är linjärt oberoende.

Homogena ekvationen har alltså endast den triviala lösningen x = 0. Vi drar slutsatsen att B^{T} B är icke-singulär och vår inhomogena ekvation B^{T} Bx = B^{T} b har unik lösning för varje högerled. Vi kallar lösningen x̂. Den ges av formeln x̂ = (B^{T} B)^{−1}B^{T} b. Vi har därmed bestämt en unik vektor P b ∈ V sådan att b − P b ⊥ V. Den kallas för den ortogonala projektionen av b på V och ges av formeln

P b = Bx̂ = B(B^{T} B)^{−1}B^{T} b.

Projektionen är en linjär operator och ges av matrisen

P = B(B^{T} B)^{−1}B^{T} .

Sats 42.10 ger en annan karaktärisering av den ortogonala projektionen: P b är det element i V som ligger närmast b.

Sats 42.10

(b − P b, v) = 0 ∀v ∈ V  ⇔  ∥b − P b∥ ≤ ∥b − v∥ ∀v ∈ V.

Läs beviset!

Exemplet: {c = [ \begin{matrix}1 \\ 1 \\ −1 \end{matrix} ] ∈ V \atop ∥b − P b∥ = ∥ [ \begin{matrix}1 − 1.286 \\ 1 − 0.571 \\ 1 − 0.857 \end{matrix} ] ∥ ≈ 0.5345 ≤ ∥b − c∥ = ∥ [ \begin{matrix}0 \\ 0 \\ 1 \end{matrix} ] ∥ = 1.}

ON-bas

Det är lättare att räkna om vi använder en ON-bas för V. Vektorerna {g_{1}, …, g_{n}} är en ON-bas om de är ortogonala och normerade dvs

(g_{i}, g_{j}) = \{ \begin{matrix}0, & i ≠ j \\ ∥g_{i}∥^{2} = 1, & i = j \end{matrix}

Med Kroneckers delta: δ_{ij} = \{ \begin{matrix}0, & i ≠ j \\ 1, & i = j \end{matrix} kan vi skriva detta (g_{i}, g_{j}) = δ_{ij}.

Vi uttrycker en allmän vektor v ∈ V med den nya basen, v = ∑_{i=1}^{n} y_{i}g_{i}. Koefficien-

terna ges av

(v, g_{j}) = (∑_{i} y_{i}g_{i}, g_{j}) = ∑_{i} y_{i} \underbrace{(g_{i}, g_{j})}_{=δ_{ij}} = y_{j}.

Vi får y_{j} = (v, g_{j}), enkelt eller hur!

Normen av vektorn uttrycks också enkelt i ON-basen. Vi har Pythagoras sats:

∥v∥^{2} = (v, v) = (∑_{i} y_{i}g_{i}, ∑_{j} y_{j}g_{j}) = ∑_{i} ∑_{j} y_{i}y_{j}(g_{i}, g_{j}) = ∑_{i} y_{i}^{2}.

Dessa räkningar kan också skrivas på matrisform. Vi bildar matrisen Q = [g_{1}, …, g_{n}], m × n. Att detta är en ON-bas betyder Q^{T} Q = I, n × n. Vektorn v ∈ V uttrycks med ON-basen, v = Qy. Koefficienterna y bestäms genom att multiplicera med Q^{T} från vänster: Q^{T} v = Q^{T} Qy = Iy = y. Pythagoras sats blir ∥v∥^{2} = v^{T} v = (Qy)^{T} (Qy) = y^{T} Q^{T} Qy = y^{T} Iy = y^{T} y = ∥y∥^{2}.

Hur skapar man en ON-bas? En ”vanlig” bas kan ortonormaliseras med hjälp av Gram-Schmidts metod. Vi visar hur det går till i vårt exempel.

Exemplet:

1. g_{1} = \frac{a_{1}}{∥a_{1}∥} = \frac{1}{\sqrt{117}} [ \begin{matrix}1 \\ 4 \\ −10 \end{matrix} ] ≈ [ \begin{matrix}0.0925 \\ 0.3698 \\ −0.9245 \end{matrix} ] (normering av första basvektorn)

V_{1} = S(a_{1}) = S(g_{1})

2. Ansats: g̃_{2} = a_{2} − α_{1}g_{1} (linjär kombination av gamla basvektorn a_{2} och den nya g_{1}).

Bestäm α_{1} så att g̃_{2} ⊥ V_{1} dvs

(a_{2} − α_{1}g_{1}, g_{1}) = 0

(a_{2}, g_{1}) − α_{1} \underbrace{(g_{1}, g_{1})}_{=1} = 0

α_{1} = (a_{2}, g_{1}) = g_{1}^{T} a_{2} = \frac{1}{\sqrt{117}}[1 4 − 10] [ \begin{matrix}2 \\ 3 \\ −5 \end{matrix} ] = \frac{64}{\sqrt{117}}

g̃_{2} = [ \begin{matrix}2 \\ 3 \\ −5 \end{matrix} ] − \frac{64}{117} [ \begin{matrix}1 \\ 4 \\ −10 \end{matrix} ] ≈ [ \begin{matrix}1.45 \\ 0.81 \\ 0.47 \end{matrix} ]

g_{2} = \frac{g̃_{2}}{∥g̃_{2}∥} = [ \begin{matrix}0.8401 \\ 0.4695 \\ 0.2718 \end{matrix} ]

Vi har nu en ON-bas Q = [g_{1}, g_{2}] = [\begin{matrix}0.0925 & 0.8401 \\ 0.3698 & 0.4695 \\ −0.9245 & 0.2718 \end{matrix}], som spänner upp samma

rum som den gamla basen B = [a_{1}, a_{2}] = [ \begin{matrix}1 & 2 \\ 4 & 3 \\ −10 & −5 \end{matrix} ].

Ekvationssystemet för ortogonal projektion blir nu: Q^{T} Qŷ = Q^{T} b, dvs ŷ = Q^{T} b, och projektionen blir P b = Qŷ = QQ^{T} b, dvs P = QQ^{T}.

Vi uttrycker vektorn c ∈ V i ON-basen: c = α_{1}g_{1} + α_{2}g_{2}, koefficienterna blir

α_{1} = (c, g_{1}) = g_{1}^{T} c ≈ [0.0925 0.3698 −0.9245] [ \begin{matrix}1 \\ 1 \\ −1 \end{matrix} ] = 1.3868,

α_{2} = (c, g_{2}) = g_{2}^{T} c ≈ [0.8401 0.4695 0.2718] [ \begin{matrix}1 \\ 1 \\ −1 \end{matrix} ] = 1.0377.

Pythagoras sats: ∥c∥^{2} = α_{1}^{2} + α_{2}^{2} ≈ 3.

Eller i matrisform:

c = Qα, α = Q^{T} c = [\begin{matrix}0.0925 & 0.3698 & −0.9245 \\ 0.8401 & 0.4695 & 0.2718 \end{matrix}] [\begin{matrix}1 \\ 1 \\ −1 \end{matrix}] ≈ [\begin{matrix}1.3868 \\ 1.0377 \end{matrix}] .

The least squares method

We now return to the original example

A = [a_{1}, a_{2}, a_{3}, a_{4}] = [ \begin{matrix}1 & 2 & 3 & 4 \\ 4 & 3 & 2 & 1 \\ −10 & −5 & 0 & 5 \end{matrix} ] , b = [ \begin{matrix}1 \\ 1 \\ 1 \end{matrix} ] .

The equation Ax = b has no solution. It is an underdetermined system because it has 4 unknowns but only 3 equations. We might also consider overdetermined systems, with more equations than unknowns.

For a system of equations Ax = b, with A of type m × n, which has no solution we might instead try to find x such that the norm of the residual b − Ax is minimal, i.e., such that

∥b − Ax∥

is minimal. From Theorem 42.10 we know that there is a unique vector P b ∈ V = R(A) which minimizes the distance between b and V. Also, since P b belongs to R(A) we can express it as at linear combination of the columns of A, i.e., P b = Ax for some x ∈ ℝ^{n}. How do we compute x? The condition that b − Ax ⊥ V leads to

\begin{matrix}(Ax, a_{j}) = (b, a_{j}), j = 1, …, n, \\ a_{j}^{T} Ax = a_{j}^{T} b, j = 1, …, n, \end{matrix}

which is the same as

A^{T} Ax = A^{T} b.

This is a symmetric system of equations with n equations and n unknowns. We have already seen that it has a solution x ∈ ℝ^{n}, but it may be non-unique. It is unique if the columns of A are linearly independent (i.e., a basis), remember that we proved this for B^{T} B before. If there are more than one solution it is customary to choose the one which is closest to the origin, i.e., the one with minimal norm ∥x∥. This is called the least-squares solution, x̂.

In Matlab: x=b\A computes the unique solution by the Gauss elimination method when A is non-singular, otherwise it automatically uses the least-squares method. If the least-squares method has multiple solutions, then it selects one with as many zero components as possible. This is ususally not the same as the one with minimal norm. To obtain this we use the so-called pseudo-inverse: x=pinv(A)*b.

In our example: Matlab calculations give

A^{T} A = [\begin{matrix}117 & 64 & 11 & −42 \\ 64 & 38 & 12 & −14 \\ 11 & 12 & 13 & 14 \\ −42 & −14 & 14 & 42 \end{matrix}] , A^{T} b = [\begin{matrix}−5 \\ 0 \\ 5 \\ 10 \end{matrix}] .

The command rref( [A'*A, A'*b] ) gives the echelon matrix

[ \begin{matrix}1 & 0 & −1 & −2 \\ 0 & 1 & 2 & 3 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \end{matrix} | \begin{matrix}−0.5429 \\ 0.9143 \\ 0 \\ 0 \end{matrix} ] ,

and we see that the system A^{T} Ax = A^{T} b has multiple solutions

x = s [ \begin{matrix}1 \\ −2 \\ 1 \\ 0 \end{matrix} ] + t [ \begin{matrix}2 \\ −3 \\ 0 \\ 1 \end{matrix} ] + [ \begin{matrix}−0.5429 \\ 0.9143 \\ 0 \\ 0 \end{matrix} ]

The command x=b\A gives

x = [\begin{matrix}0.0667 \\ 0 \\ 0 \\ 0.3048 \end{matrix}] .

Here we have s = 0, t = 0.3048. The command x=pinv(A)*b gives

x̂ = [\begin{matrix}−0.0143 \\ 0.0571 \\ 0.1286 \\ 0.2000 \end{matrix}] .

Here we have s = 0.1286, t = 0.2000.

2004-02-07 /stig