In [1]:
#The Hilbert matrix!
#=
Suppose that we are interested in representing a function
f(x) over the interval $[0,1]$ with a polynomial approximation.
This problem was considered by Hilbert around 1895.
https://blogs.mathworks.com/cleve/2017/06/07/hilbert-matrices-2/
This is a standard computation in numerical analysis where
we seek some way to represent $f(x)$ on the computer, and
polynomials can represent any continuous function to arbitrary
accuracy. (Internally, this is often closely related to how
the computer computes things likes $\sin(x)$)
Back to the problem at hand, if we approximate $f(x)$
with $p(x; c) = \sum_{j=1}^c c_i x^{j-1}$
in the sense that we want:
$\int_0^1 (p(x; c) - f(x))^2 dx$ to be as small as possible,
then we can show that is equivalent to minimizing
$\int_0^1 (p(x; c) - f(x))^2 dx = c^T A c - 2 c^T b + \theta$.
where the term $A_{ij} = \int_0^1 x^{i-1} x^{j-1} \, dx$,
$b_i = \int_0^1 f(x) x^{i-1} \, dx$, and
$\theta = \int_0^t f(x)^2 \, dx$.
Just to be clear, here we have that
$c^T A c = \int_0^1 p(x; c)^2\, dx > 0$ as long as $c \not= 0$.
hence the matrix $A$ is positive definite!
Thus, the solution to
$\min c^T A c - 2 c^T b + \theta$ is given by the linear system
$ A c = b$.
The entries $A_{ij} = \int_0^1 x^{i-1} x^{j-1} \, dx = 1/(i+j-1)$.
If we are interested in approximating the function $f(x) = \sqrt(x)$,
then $b_j = \int_0^1 x^{1/2} x^{j-1} \, dx = 1/(j+1/2)$
=#
# Let's solve this!
function hilbert(n::Integer)
A = zeros(n,n)
for j=1:n
for i=1:n
A[i,j] = 1/(i+j-1)
end
end
return A
end
hilbert(5)
Out[1]:
In [2]:
##
function approxprob(k)
A = hilbert(k)
b = zeros(k)
for j=1:k
b[j] = 1/(j+1/2)
end
return A,b
end
A, b = approxprob(10)
c = A\b
Out[2]:
In [3]:
##
using Polynomials
using Plots
pyplot()
p = Poly(c)
plot(x -> p(x), range(0,stop=1, length=50))
gui()
In [4]:
##
using LinearAlgebra
plot()
for k in [2,5,10,20,50,100,200]
A,b = approxprob(k)
c = A\b
p = Poly(c)
@show norm(A*c - b)
plot!(x -> p(x), range(0,stop=1, length=50), lab="$k")
end
gui()
In [5]:
## This appears to work great!
In [6]:
## Let's make a small perturbation
using LinearAlgebra
plot()
for k in [2,5,10,20,50,100,200]
A,b = approxprob(k)
b .+= 1e-11*randn(size(A,1))
c = A\b
p = Poly(c)
@show norm(A*c - b)
plot!(x -> p(x), range(0,stop=1, length=50), lab="$k")
end
gui()
In [7]:
##
In [8]:
##
# But if you try solving with this system with a random right hand side...
for k in [2,5,10,20,50,100,200]
A,b = approxprob(k)
b = randn(k)
c = A\b
@show k, norm(A*c-b)/norm(b)
end
# it doesn't work at all!
# so, here's a question, is there a simple to evaluate
# function f(x) (e.g. int_0^1 f(x) x^j dx has a
# simple expression, such that this shows the ill-conditioning?)
In [9]:
## Added in class
for k in [2,5,10,20,50,100,200]
A,b = approxprob(k)
b = randn(k)
Q,R = qr(A)
c = R\(Q'*b)
@show k, norm(A*c-b)/norm(b)
end