MATH 417 Lecture 25

From Notes
Jump to navigation Jump to search

« previous | Tuesday, April 22, 2014 | next »


Section 8.1 and 8.2

Least Squares approximation

f(x)≈a1ϕ1(x)+a2ϕ2(x)+…+anϕn(x)

usually n=2 (linear) or n=3 (quadratic)

⟨f,g⟩ is given, so we wish to minimize ‖f−p‖2=⟨f−p,f−p⟩

In 8.1, ⟨f,g⟩=∑i=1nwifigi, where wi>0 is a probability weight distribution.

In 8.2, ⟨f,g⟩=∫abw(x)f(x)g(x)

Common definitions of w(x):

  • w(x)=1 (classic definition)
  • w(x)=11+x2 on [−∞,∞]
  • w(x)=1−x2 on [−1,1]
  • w(x)=11−x2

Basic Concept

We compute Ψ(a1,a2,…,an)=∫01|f(x)−p(x)|2w(x)dx and find the minimum. To find this minimum, we solve for ∇Ψ=0. This gives a system of normal equations:

A=[⟨ϕ1,ϕ1⟩⟨ϕ1,ϕ2⟩…⟨ϕ1,ϕn⟩⟨ϕ2,ϕ1⟩⟨ϕ2,ϕ2⟩…⟨ϕ2,ϕn⟩⋮⋮⋱⋮⟨ϕn,ϕ1⟩⟨ϕ2,ϕn⟩…⟨ϕn,ϕn⟩]

This matrix is positive definite, so it has an inverse.

b=[⟨ϕ1,f⟩⟨ϕ2,f⟩⋮⟨ϕn,f⟩]

This system is solvable, but it's a pain to solve.

Original function f(x)=x3.
Best Quadratic fit p(x)=120−35x+32x2.

For example, find the best fit to x3 by p(x)=a1+a2x+a3x2 on [0,1]

[11213121314131415][a1a2a3]=[⟨1,x3⟩⟨x,x3⟩⟨x2,x3⟩]


Simplifying Normal Equations

What if A were a diagonal matrix? That is, what if ⟨ϕi,ϕj⟩=0 for i≠j?

Then ϕ1, ϕ2, etc. would be an orthogonal basis for our regression space.


Continuing with our example, let p=c1ϕ1+c2ϕ2+c3ϕ(3), where ϕ1(x)=1, ϕ2(x)=2x−1, and ϕ3(x)=3(2x−1)2−1. {ϕ1,ϕ2,ϕ3} is an orthogonoal basis for our regression space.

[⟨1,1⟩000⟨2x−1,2x−1⟩000⟨3(2x−1)2−1,3(2x−1)2−1⟩][c1c2c3]=[⟨1,x3⟩⟨2x−1,x3⟩⟨3(2x−1)2−1,x3⟩]


This gives a much easier system to solve:

[10001300045][c1c2c3]=[1425−14110]


There are two ways to find an orthogonal basis for p

  1. Legendre polynomials
  2. Gram-Schmidt

Legendre Polynomials

Ln(x)=dndxn[(1−x2)n]

For n∈{0,1,2}, we have:

  • L0(x)=1,
  • L1(x)=x, and
  • L2(x)=3x2−1

These only work on [−1,1], so we shift and scale to [0,1] by substituting y=x+12 (or x=2y−1). Performing this substitution gives our ϕ's above.


Gram-Schmidt

This will always give an orthonormal basis, so ⟨ϕi,ϕj⟩=δij:

def ip(f, g, a, b):
    return integrate(f*g, x, a, b)

def norm(f, a, b):
    return sqrt(ip(f, f, a, b))

def normalize(f, a, b):
    return f/norm(f, a, b)

def gram_schmidt(basis, a, b):
    on_basis = [normalize(basis[0], a, b)]
    for i in xrange(1, len(basis)):
        p = sum(ip(basis[i], e, a, b) * e for e in on_basis)
        next_e = normalize(basis[i] - p, a, b)
        on_basis.append(next_e)
    return map(simplify, on_basis)
  • φ1(x)=1
  • φ2(x)=3(2x−1)
  • φ3(x)=5(6x2−6x+1)


Now we have

[100010001][α1α2α3]=[143320520]


Quiz Discussion: Matrix Norms

‖A‖p=max‖x‖p=1‖Ax‖p

This is horrible. Bojan Popov


Find ‖A‖1, ‖A‖2, and ‖A‖∞ for

A=[10−2020−103]

‖A‖∞ is the max row sum, so ‖A‖∞=4

‖A‖1 is just the max column sum, so ‖A‖1=5

‖A‖2=ρ(ATA), so

  1. ATA=[20−5040−5013],
  2. ρ(ATA)=maxiλi=max⁡{4,15+2212,15−2212}=15+2212
  3. ρ(ATA)=15+2212