Linear Least Squares Fitting

Introduction
Linear Regression
Linear Approximation Theory: an introduction

Introduction

We are likely all familiar with the concept of curve-fitting — sometimes, we, as humans, can recognize certain patterns in data about which the data itself is agnostic. If we look at a chart like the following: We think immediately of a linear relationship. But what tools can we use to precisely determine a line from these individual points of data? Intuitively speaking, we can all figure this out pretty easily — we want to find the line that is minimally far from every single one of these points. The way this is actually implemented is through a method known as least squares regression. This form of regression is very powerful, and is widely used in applications from signals processing to economics research. Least squares fitting finds the best curve to fit a set of points through minimizing the sum of the squares of the offsets of each point from the curve.

Linear Regression

Suppose that we have measured three data points $P = {(- 1, 3), (1, 1), (2, 2)}$ , and that our model for these data asserts that the points should lie on a line. Of course, these three points do not actually lie on a single line, but this could be due to errors in our measurement. How do we predict which line they are supposed to lie on?

The general equation for a (non-vertical) line is $y = a_{0} + a_{1} x$ . If our three data points were to lie on this line, then the following equations would be satisfied: ${\begin{cases} a_{0} - a_{1} = 3 \\ a_{0} + a_{1} = 1 \\ a_{0} + 2 a_{1} = 2 \end{cases}$ In order to find the best-fit line, we try to solve the above equations in the unknowns $a_{0}$ and $a_{1}$ . Putting our linear equations into matrix form, we are trying to solve $A x = b$ where $A = (\begin{array}{cc} 1 & - 1 \\ 1 & 1 \\ 1 & 2 \end{array}) x = (\begin{matrix} a_{0} \\ a_{1} \end{matrix}) b = (\begin{matrix} 3 \\ 1 \\ 2 \end{matrix})$ As the three points do not actually lie on a line, there is no actual solution. So instead we have to compute an approximate solution.

Consequently, we have to answer to the following important question:
Suppose that $A x = b$ does not have a solution. What is the best approximate solution?

For our purposes, the best approximate solution is called the least-squares solution. Least squares fitting finds the best curve to fit a set of points through minimizing the sum of the squares of the offsets of each point from the curve. We use the squares of the offsets to ensure the differentiability property, which allows us to linear tools from calculus in analyzing the data. We will present two methods for finding least-squares solutions, and we will give several applications to best-fit problems. We begin by clarifying exactly what we will mean by a “best approximate solution” to an inconsistent matrix equation $A x = b$ .

Let

A

be a

m \times n

(with

m > n

) matrix of

M_{m, n} (R)

and

b

be a vector of

R^{m}

. A least-squares solution of the matrix equation

A x = b

is a vector

x^{*}

R^{n}

such that

for all x \in R^{n} {‖ {A . x}^{*} - b ‖}_{2} \leq {‖ A . x - b ‖}_{2}

The term “least squares” comes from the fact that ${‖ A . x - b ‖}_{2}$ is the square root of the sum of the squares of the entries of the vector $A . x - b$ . So a least-squares solution minimizes the sum of the squares of the differences between the entries of $A x$ and $b$ . In other words, a least-squares solution solves the equation $A x = b$ as closely as possible, in the sense that the sum of the squares of the difference $A . x - b$ is minimized.

Let

A

be a

m \times n

(with

m > n

) matrix of

M_{m, n} (R)

and

b

be a vector of

R^{m}

. The least-squares solutions of

A x = b

are the solutions of the matrix equation

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

where

A^{T}

is the transpose matrix of

A

. Moreover, if

rank (A) = n

then

A x = b

has a unique least-squares solution

x^{*}

with

x^{*} = A^{+} . b

where

A^{+} = (A^{T} . A)^{- 1} . A^{T}

is called the pseudo inverse of matrix

A

Example: Let $y = a_{0} + a_{1} x$ be a linear model of the data $P = {(0, 1), (2, 1.9), (4, 3.2)}$ . Then we have to solve the overdetermined linear system $A x = b$ by using a least squares solution where $A^{T} = [\begin{array}{ccc} 1 & 1 & 1 \\ 0 & 2 & 4 \end{array}], b^{T} = [1, 1.9, 3.2]$ As the rank of the matrix is equal to $2$ , by using the previous Theorem, it yields: $\hat{x} = A^{+} b = (A^{T} A)^{- 1} A^{T} b = [0.93, 0.55]^{T}$ Consequently, $y = 0.93 + 0.55 x$ is the unique least squares solution to this problem.

The least squares method is based on the residual, the difference between the predicted value by the model

A x

and the actual value

b

given in the dataset:

r = A x - b .

To find the optimal function parameters in

x

, we have to minimize the squared error

F (x) = | | r | |_{2}^{2} = r^{T} r = (A x - b)^{T} (A x - b)

F (x) = x^{T} A^{T} . A x - 2 x^{T} A^{T} b + b^{T} b

A necessary condition is to minimized this quadratic function:

F^{'} (x) = 2 (x^{T} A^{T} A - b^{T} A) = 0

rank (A) = n

then matrix

A^{T} A

is positive definite i.e. this matrix is invertible:

for any vector x \in R^{n} x^{T} A^{T} A x = {‖ A x ‖}_{2} \geq 0

rank (A) = n

it yields that

x^{T} A^{T} A x = 0

implies that

A x = 0

and

x = 0

.
Consequently,

A^{T} A

is invertible and

x^{*} = (A^{T} A)^{- 1} A^{T} b = A^{+} b

Sufficient condition:
if

x^{*}

satifies

{A^{T} . A . x}^{*} = A^{T} . b

then for any

x \in R^{n}

we have

F (x) - F (x^{*}) = (A x - b)^{T} (A x - b) - ({A x}^{*} - b)^{T} ({A x}^{*} - b),

F (x) - F (x^{*}) = (x - x^{*})^{T} A^{T} A (x - x^{*})

After some calculations, it yields

(x - x^{*})^{T} A^{T} A (x - x^{*}) = {‖ A (x - x^{*}) ‖}_{2}^{2} \geq 0

Consequently, for any

x \in R^{n}

we have

F (x) \geq F (x^{*})

. Moreover if

rank (A) = n

then

F (x) = F (x^{*})

implies that

{‖ A (x - x^{*}) ‖}_{2} = 0

and

x = x^{*}

. Hence,

x^{*}

is the unique minimum of

F

You may have noticed that we have been careful to say “the least-squares solutions” in the plural, and “a least-squares solution” using the indefinite article. This is because a least-squares solution need not be unique: indeed, if the columns of $A$ are linearly dependent, then the least-squares solution of $A x = b$ has infinitely many solutions. The previous theorem, gives us a criteria for uniqueness.

Exercice [PY-1]

We wish to calculate the least-square approximation of the previous data $P = {(- 1, 3), (1, 1), (2, 2)}$ with a polynomial function of degree $n = 1$ .
Firstly, write this problem as an interpolation problem of the form $A . x = b$ . We then obtain an overdetermined system. We propose to solve it by minimizing the following function: $min_{(a_{0}, a_{1}) \in R^{2}} F (a_{0}, a_{1}) = | | A . x - b | |^{2}$ Show that $| | A . x - b | |^{2} = \sum_{j = 1}^{3} {(a_{0} + a_{1} x_{j} - y_{j})}^{2}$
By writting all the partial derivatives of $F$ and by solving $\frac{\partial F}{\partial a_{0}} = 0$ and $\frac{\partial F}{\partial a_{1}} = 0$ show that we have to solve $A^{T} A . x = A^{T} b$ .
Draw the data points and the linear function solution to this system.

f (x) = \cos (x)

[0, 2 π]

p = 10

p = 20

P = {(x_{j} = \frac{(j - 1)}{p} π, y_{j} = f (x_{j})) f o r j = 1, . . ., p}

Calculate the best least-squares polynomial approximation $P_{4} (x) = a_{0} + a_{1} x + a_{2} x^{2} + a_{3} x^{3} + a_{4} x^{4}$ of degree equal or less than $n = 4$ of these data points $P$ .
Print the maximum error.
Draw the data points and the polynomial solution $P_{4} (x)$ .
In this part we want to minimize the sum of the squares of the relative errors: $ε_{j} = \frac{P_{4} (x_{j}) - y_{j}}{y_{j}}$ For this we want to minimize the following function: $\begin{aligned} F (a_{0}, a_{1}, \dots, a_{4}) & = \sum_{j = 1}^{p} ε_{j}^{2} \\ = \sum_{j = 1}^{p} w_{j} {(P_{4} (x_{j}) - y_{j})}^{2} \end{aligned}$ where the values $w_{j} = \frac{1}{{(y_{j})}^{2}}$ are called the weights of the data points. We assume that $y_{j} \neq 0$ . Consequently, the data with high values $y_{j}$ are less sensitive in the minimization process.
Show that, in this case, we have to solve the following linear system: $A^{T} . W . A . x = A^{T} . W . b$ where $W$ is a $p \times p$ diagonal matrix with the $w_{j}$ as diagonal values.

Firstly, we calculate the partial derivatives:

\frac{\partial F}{\partial a_{0}} = 2 \times \sum_{j = 1}^{3} 1 \times (a_{0} + a_{1} x_{j} - y_{j}) = 0

and

\frac{\partial F}{\partial a_{1}} = 2 \times \sum_{j = 1}^{3} x_{j} (a_{0} + a_{1} x_{j} - y_{j}) = 0

We then obtain the following linear system:

{\begin{cases} a_{0} \sum_{j = 1}^{3} 1 + a_{1} \sum_{j = 1}^{3} x_{j} & = \sum_{j = 1}^{3} y_{j} \\ a_{0} \sum_{j = 1}^{3} x_{j} + a_{1} \sum_{j = 1}^{3} x_{j}^{2} & = \sum_{j = 1}^{3} x_{j} y_{j} \end{cases}

Consequently we have to solve

B . x = c

where

B = [\begin{array}{ll} 3 & \sum_{j = 1}^{3} x_{j} \\ \sum_{j = 1}^{3} x_{j} & \sum_{j = 1}^{3} x_{j}^{2} \end{array}], c = [\begin{array}{l} \sum_{j = 1}^{3} y_{j} \\ \sum_{j = 1}^{3} x_{j} . y_{j} \end{array}]

and

x^{T} = [a_{0}, a_{1}]

.
By direct calculations we can show that

B = A^{T} A

and

c = A^{T} b

For

n = 4

we have to solve:

{\begin{cases} a_{0} . p + a_{1} \sum_{j = 1}^{p} x_{j} + a_{2} \sum_{j = 1}^{p} x_{j}^{2} + a_{3} \sum_{j = 1}^{p} x_{j}^{3} + a_{4} \sum_{j = 1}^{p} x_{j}^{4} & = \sum_{j = 1}^{p} y_{j} \\ a_{0} \sum_{j = 1}^{p} x_{j} + a_{1} \sum_{j = 1}^{p} x_{j}^{2} + a_{2} \sum_{j = 1}^{p} x_{j}^{3} + a_{3} \sum_{j = 1}^{p} x_{j}^{4} + a_{4} \sum_{j = 1}^{p} x_{j}^{5} & = \sum_{j = 1}^{p} x_{j} y_{j} \\ a_{0} \sum_{j = 1}^{p} x_{j}^{2} + a_{1} \sum_{j = 1}^{p} x_{j}^{3} + a_{2} \sum_{j = 1}^{p} x_{j}^{4} + a_{3} \sum_{j = 1}^{p} x_{j}^{5} + a_{4} \sum_{j = 1}^{p} x_{j}^{6} & = \sum_{j = 1}^{p} x_{j}^{2} y_{j} \\ a_{0} \sum_{j = 1}^{p} x_{j}^{3} + a_{1} \sum_{j = 1}^{p} x_{j}^{4} + a_{2} \sum_{j = 1}^{p} x_{j}^{5} + a_{3} \sum_{j = 1}^{p} x_{j}^{6} + a_{4} \sum_{j = 1}^{p} x_{j}^{7} & = \sum_{j = 1}^{p} x_{j}^{3} y_{j} \\ a_{0} \sum_{j = 1}^{p} x_{j}^{4} + a_{1} \sum_{j = 1}^{p} x_{j}^{5} + a_{2} \sum_{j = 1}^{p} x_{j}^{6} + a_{3} \sum_{j = 1}^{p} x_{j}^{7} + a_{4} \sum_{j = 1}^{p} x_{j}^{8} & = \sum_{j = 1}^{p} x_{j}^{4} y_{j} \end{cases}

F (a_{0}, a_{1}, \dots, a_{4}) = \sum_{j = 1}^{p} w_{j} {(P_{4} (x_{j}) - y_{j})}^{2}

then for any

k \in {0, 1, 2, 3, 4}

\frac{\partial F}{\partial a_{k}} = 2 \times \sum_{j = 1}^{p} x_{j}^{k} w_{j} \times (a_{0} + a_{1} x_{j} + . . . + a_{4} x_{j}^{4} - y_{j}) = 0

It yields the following equations for any

k \in {0, 1, 2, 3, 4}

a_{0} \sum_{j = 1}^{p} w_{j} x_{j}^{k} + a_{1} \sum_{j = 1}^{p} w_{j} x_{j}^{k + 1} + . . . + a_{4} \sum_{j = 1}^{p} w_{j} x_{j}^{k + 4} = \sum_{j = 1}^{p} w_{j} x_{j}^{k} y_{j}

which can be written:

A^{T} . W . A . x = A^{T} . W . b

with

W = diag (w_{1}, w_{2}, . . ., w_{p})

import numpy as np
import matplotlib.pyplot as plt
import numpy.linalg as LA

print("*******************") 
print("PY-1")
print("*******************") 



# création de la liste des données à interpoler
def f(x):
    return np.abs(x)+1

# Horner Algorithm
def Horner(x,coeff):
    p=len(coeff)
    s=coeff[p-1]
    for i in range(p-1):
        s=x*s+coeff[p-2-i]
    return s

# nombre de points
p=7
# degré du polynôme
n=4
# intervalle
a=-np.pi
b= np.pi

X=np.linspace(a,b,p)
Y=f(X)

# Création de la matrice d'interpolation
A=np.zeros((p,n+1))

for i in range(p):
    for j in range(n+1):
        A[i,j]=X[i]**j
        
#print("matrice interpolation: ")
#print(np.round(A,1))
# résolution
coeff=LA.solve(np.dot(np.transpose(A),A),np.dot(np.transpose(A),Y))
#print("mat= ",np.dot(np.transpose(A),A)," b= ",np.dot(np.transpose(A),Y))
print("solution= ",coeff)

err1 = [np.abs(Horner(X[i],coeff)-Y[i]) for i in range(p)]
print("\nmax error data points = ",np.max(err1))

t=np.linspace(a,b,200)
err2 = [np.abs(Horner(t[i],coeff)-f(t[i])) for i in range(len(t))]
print("\nmax error data |polynom -function| = ",np.max(err2))

print("\nRelative error minimization")

w = [1/Y[i]**2 for i in range(p)]
W = np.diag(w)

# résolution
coeffr=LA.solve(np.dot(np.transpose(A),np.dot(W,A)),np.dot(np.transpose(A),np.dot(W,Y)))
print("solution= ",coeffr)
err1 = [np.abs(Horner(X[i],coeffr)-Y[i]) for i in range(p)]
print("\nmax error data points = ",np.max(err1))

t=np.linspace(a,b,200)
err2 = [np.abs(Horner(t[i],coeffr)-f(t[i])) for i in range(len(t))]
print("\nmax error data |polynom -function| = ",np.max(err2))

plt.clf()
p1=plt.plot(X,Y,'bo')
p2=plt.plot(t,f(t),'k',label='function')
p3=plt.plot(t,Horner(t,x),'r',label="l.s. polynom")
p3=plt.plot(t,Horner(t,xr),'g',label="l.s.r polynom")
plt.legend(loc='lower right')
plt.show()
plt.savefig('PY1.png')

Exercice [PY-2]

Gauss invented the method of least squares to find a best-fit ellipse: he correctly predicted the (elliptical) orbit of the asteroid Ceres as it passed behind the sun in 1801.

Find the best least square ellipse through the points $P = {(0, 2), (2, 1), (1, - 1), (- 1, - 2), (- 3, 1), (- 1, - 1)}$ . What quantity is being minimized?
Solve this problem with the numpy python library and draw the best least-squares ellipse solution. For this, you can use the contour function of matplotlib.pyplot

The general equation for an ellipse (actually, for a nondegenerate conic section) is

f (x, y) = x^{2} + B y^{2} + C x y + D x + E y + F = 0.

This is an implicit equation: the ellipse is the set of all solutions of the equation. To say that our data points lie on the ellipse means that the above equation is satisfied for the given values of

x

and

y

. To put this in matrix form, we move the constant terms to the right-hand side of the equals sign; then we can write this as

A x = b

and solve it by using the normal equations.

print("\n***************")
print("Exercice PY-2")
print("***************\n")

def ellipse(x,y,c):
    res = x*x+c[0]*y*y+c[1]*x*y+c[2]*x+c[3]*y+c[4]
    return res

P=[[0,2],[2,1],[1,-1],[-1,-2],[-3,1],[-1,1]]
xp = np.transpose(P)[0]
yp = np.transpose(P)[1]  

# Approximation au sens des moindres carrées de la conique
A=[]
b=[]
for i in range(len(P)):
    x=P[i][0]
    y=P[i][1]
    A.append([y*y,x*y,x,y,1])
    b.append(-x*x)
coeff=LA.solve(np.dot(np.transpose(A),A),np.dot(np.transpose(A),b))

#calcul de la grille des valeurs et de la conique
delta = 0.025
x = np.arange(-5.0, 5.0, delta)
y = np.arange(-3.0, 3.0, delta)
X, Y = np.meshgrid(x, y)
Z=ellipse(X,Y,coeff)


plt.clf()
#tracé des points
p1=plt.plot(xp,yp,'bo')
CS=plt.contour(X, Y, Z,levels = np.arange(0.0, 1.0, 5.0),extend='both')
plt.clabel(CS,inline=1,fontsize=10)
plt.title("Least Squares Ellipse")
plt.savefig('ellipse.png')

Segmented Linear Regression

Segmented linear regression (SLR) addresses this issue by offering piecewise linear approximation of a given dataset [2]. It splits the dataset into a list of subsets with adjacent ranges and then for each range finds linear regression, which normally has much better accuracy than one line regression for the whole dataset. It offers an approximation of the dataset by a sequence of line segments. In this context, segmented linear regression can be viewed as a decomposition of a given large dataset into a relatively small set of simple objects that provide compact but approximate representation of the given dataset within a specified accuracy. Segmented linear regression provides quite robust approximation. This strength comes from the fact that linear regression does not impose constraints on ends of a line segment. For this reason, in segmented linear regression, two consecutive line segments normally do not have a common end point. This property is important for practical applications, since it enables us to analyse discontinuities in temporal data: sharp and sudden changes in the magnitude of a signal, switching between quasi-stable states. This method allows us to understand the dynamics of a system or an effect. For comparison, high degree polynomial regression provides a continuous approximation [3]. Its interpolation and extrapolation suffer from overfitting oscillations [4]. Polynomial regression is not a reliable predictor in applications mentioned earlier. Despite the simple idea, segmented linear regression is not a trivial algorithm when it comes to its implementation. Unlike linear regression, the number of segments in an acceptable approximation is unknown. The specification of the number of segments is not a practical option, since a blind choice produces usually a result that badly represents a given dataset. Similarly, uniform splitting a given dataset into equal-sized subsets is not effective for achieving an optimal approximation result. Instead of these options, we place the restriction of the allowed deviation and attempt to construct an approximation that guarantees the required accuracy with as few line segments as possible. The line segments detected by the algorithm reveal essential features of a given dataset. The deviation at a specified x value is defined as the absolute difference between original and approximation y values. The accuracy of linear regression in a range and the approximation error are defined as the maximum of all of the values of the deviations in the range. The maximum allowed deviation (tolerance) is an important user interface parameter of the model. It controls the total number and lengths of segments in a computed segmented linear regression. If input value of this parameter is very large, there is no splitting of a given dataset. The result is just one line segment of linear regression. If input value is very small, then the result may have many short segments. It is interesting to note that within this formulation, the problem of segmented linear regression is similar to the problem of polyline simplification. The solution to this problem provides the algorithm of Ramer, Douglas and Peucker [5]. Nevertheless, there are important reasons why polyline simplification is not the best choice for noisy data. This algorithm gives first priority to data points the most distant from already constructed approximation line segments. In the context of this algorithm, such data points are the most significant. The issue is that in the presence of noise, the significant points normally arise from high frequency fluctuations. The biggest damage is caused by outliers. For this reason, noise can have substantial effect on the accuracy of polyline simplification. In piecewise linear regression, each line segment is balanced against the effect of noise. The short term fluctuations are removed to show significant long term features and trends in original data.

Linear Approximation Theory: an introduction

Approximation theory is concerned with the ability to approximate functions by simpler and more easily calculated functions. The first question we ask in approximation theory concerns the possibility of approximation. Is the given family of functions from which we plan to approximate dense in the set of functions we wish to approximate? In this section we survey some of the main density results and density methods.

Let us consider the following function $f (x) = | x |$ , we want to find a polynom function $P_{n} (x)$ of degree $n$ which closely approximates $f$ . One way to choose this polynom is to minimized the "area" or the "distance" between $f$ and $P_{n}$ .

Let $E = C^{0} ([a, b])$ be the vector space of continuous functions on $[a, b]$ (with $a < b$ ) into $R$ . We can easily show that the polynomial vector space $P_{n}$ belongs to $C^{0} ([a, b])$ .

We denote by $⟨ f, g ⟩ : E \times E \to R$ the scalar product of any functions $f$ and $g$ belonging to $E$ . It satisfies $\begin{array}{lll} \cdot ⟨ f + h, g ⟩ = ⟨ f, g ⟩ + ⟨ h, g ⟩ & \cdot for any λ \in R ⟨ λ f, g ⟩ = λ . ⟨ f, g ⟩ \\ \cdot ⟨ f, g ⟩ = ⟨ g, f ⟩ & \cdot ⟨ f, f ⟩ \geq 0 & \cdot ⟨ f, f ⟩ = 0 \Rightarrow f = 0 \end{array}$

We say that such scalar product on

E

is a symetric bilinear form which is positive definite.

The scalar product can be defined by:

⟨ f, g ⟩ = \int_{a}^{b} w (x) f (x) g (x) d x

where

w (x)

is a "weighted" function. For example we can choose:

$w (x) = 1$ or $w (x) = \frac{1}{\sqrt{1 - x^{2}}}$ on $] - 1, 1 [$
$w (x) = e^{- x^{2}}$ on $] - \infty, \infty [$

So as to define a polynomial approximation

P_{n}

of a continuous function

f

, we need a metric or a distance between

f

and

P_{n}

. For this, we define a norm on

E

Let $E$ be the vector space of continuous functions on $[a, b]$ and $⟨ ., . ⟩$ be a scalar product on $E$ then ${‖ f ‖}_{2} = \sqrt{⟨ f, f ⟩}$ is a norm on $E$ .

In our previous example, we have:

{‖ f ‖}_{2} = {(\int_{a}^{b} w (x) {| f (x) |}^{2} d x)}^{\frac{1}{2}}

. Consequently, the distance between two functions is given by

d (f, g) = {‖ f - g ‖}_{2}

Let us proof the Cauchy -Schwarz inequality

for any f, g \in E, we have | ⟨ f, g ⟩ | \leq {‖ f ‖}_{2} {‖ g ‖}_{2}

For any

λ \in R

and by using the scalar product properties, it yields:

⟨ f + λ g, f + λ g ⟩ = ⟨ f, f ⟩ + 2 λ ⟨ f, g ⟩ + λ^{2} ⟨ g, g ⟩ \geq 0

⟨ g, g ⟩ \neq 0

, we set

λ = - \frac{⟨ f, g ⟩}{⟨ g, g ⟩}

then we have

⟨ f, f ⟩ - \frac{(⟨ f, g ⟩)^{2}}{⟨ g, g ⟩} \geq 0

and the inequality follows.
If

⟨ g, g ⟩ = 0

then for any

λ \in R

we have

⟨ f, f ⟩ + 2 λ ⟨ f, g ⟩ \geq 0

. As

⟨ f, f ⟩ \geq 0

, it implies that

⟨ f, g ⟩

must be equal to

0

. The inequality is then also satisfied.

We can then show that

{‖ f ‖}_{2}

is a norm.

${‖ f + g ‖}_{2}^{2} = ⟨ f, f ⟩ + 2 ⟨ f, g ⟩ + ⟨ g, g ⟩ = {‖ f ‖}_{2}^{2} + 2 ⟨ f, g ⟩ + {‖ g ‖}_{2}^{2}$ By using the Cauchy-Schwarz inequality it yields: ${‖ f + g ‖}_{2}^{2} \leq {‖ f ‖}_{2}^{2} + 2 {‖ f ‖}_{2} {‖ g ‖}_{2} + {‖ g ‖}_{2}^{2} = ({‖ f ‖}_{2} + {‖ g ‖}_{2})^{2}$ Consequently ${‖ f + g ‖}_{2} \leq {‖ f ‖}_{2} + {‖ g ‖}_{2} .$

The other norm properties directly follow from the scalar product properties:

{‖ f ‖}_{2} = (⟨ f, f ⟩)^{\frac{1}{2}} \geq 0

{‖ f ‖}_{2} = 0 \Rightarrow ⟨ f, f ⟩ = 0 \Rightarrow f = 0

and

{‖ λ f ‖}_{2} = (λ^{2} ⟨ f, f ⟩)^{\frac{1}{2}} = | λ | . {‖ f ‖}_{2}

Exercice [M-1]

Show that for any $f \in E$ , ${‖ f ‖}_{2} = {(\int_{a}^{b} {| f (x) |}^{2} d x)}^{\frac{1}{2}}$ is a norm.
Show that ${‖ f ‖}_{1} = \int_{a}^{b} | f (x) | d x$ is also a norm. Show that we cannot associate a scalar product to this norm.
We set $a = 0$ and $b = 1$ . We consider the following functions $f (x) = x^{2}$ and $P_{λ} (x) = λ x$ where $λ \in [0, 1]$ . Calculate the optimal values of $λ$ in the following minimization problems: $min_{λ \in [0, 1]} {‖ f - P_{λ} ‖}_{1} and min_{λ \in [0, 1]} {‖ f - P_{λ} ‖}_{2}$

We start with the

1

-norm demonstration:

- We have ${‖ f ‖}_{1} = \int_{a}^{b} | f (x) | d x \geq 0$ as the integration of a positive function is positive. The fact that $a < b$ is very important here.
- Moreover we have that ${‖ 0 ‖}_{1} = 0$ and ${‖ λ f ‖}_{1} = | λ | . \int_{a}^{b} | f (x) | d x = | λ | . {‖ f ‖}_{1}$ due to the integral properties.
- The triangular inequality is straightforward by considering for any $x \in [a, b]$ the following inequality: $| f (x) + g (x) | \leq | f (x) | + | g (x) |$ Consequently ${‖ f + g ‖}_{1} \leq \int_{a}^{b} (| f (x) | + | g (x) |) d x = \int_{a}^{b} | f (x) | d x + \int_{a}^{b} | g (x) | d x$ by the integral monotony property. It yields ${‖ f + g ‖}_{1} \leq {‖ f ‖}_{1} + {‖ g ‖}_{1}$ .
- To conclude, we have to show that for any positive and continuous function $g$ on $[a, b]$ that $\int_{a}^{b} g (x) d x = 0 \Rightarrow for any x \in [a, b] g (x) = 0.$ We will proof it by contrapositive (we assume $n o n (Q)$ and show that this implies $n o n (P)$ ). We set $g (x) = | f (x) |$ . We assume that there exists $x_{0} \in [a, b]$ such that $g (x_{0}) = γ > 0$ . As function $g$ is continuous, by setting $ϵ = \frac{γ}{2} > 0$ we know that: $there exists η > 0 such that for any x \in [x_{0} - η, x_{0} + η], | g (x) - γ | \leq \frac{γ}{2} .$ Consequently for any $x \in [x_{0} - η, x_{0} + η]$ , $\frac{γ}{2} \leq g (x) \leq \frac{3 γ}{2}$ . As $g \geq 0$ on $[a, b]$ it yields $\int_{a}^{b} g (x) d x \geq \int_{x_{0} - η}^{x_{0} + η} g (x) d x \geq \int_{x_{0} - η}^{x_{0} + η} \frac{γ}{2} d x = γ . η > 0$ . This achieves the proof.
Remark. We can adapt the interval to $[a, a + η]$ (resp. $[b - η, b]$ ) if $x_{0} = a$ (resp. $x_{0} = b$ ).

The

2

-norm case:

It is easy to show that ${‖ f ‖}_{2} = (ϕ (f, f))^{\frac{1}{2}}$ where $ϕ (f, f) = \int_{a}^{b} {| f (x) |}^{2} d x$ . If $ϕ (f, f)$ is a scalar product it then follows directly by the previous Theorem that ${‖ f ‖}_{2}$ is a norm. Let us show the more difficult property: $ϕ (f, f) = 0$ implies that $f = 0$ . For this, we set $g (x) = | f (x) |^{2}$ . As $ϕ (f, f) = \int_{a}^{b} g (x) d x$ then by using the previous proof, we obtain that $ϕ (f, f) = 0$ implies that $g (x) = | f (x) |^{2} = 0$ and consequently $f = 0$ .

We set

ϕ (x) = | f (x) - λ . x | = x | x - λ |

. We study this function and we find that

λ_{1} = \frac{1}{\sqrt{2}}

and

λ_{2} = \frac{3}{4}

Let $P_{n}$ be a polynom of $P_{n}$ . The polynom $P_{n}$ is the best least square approximation of $f \in [a, b]$ if ${‖ f - P_{n} ‖}_{2} = min_{Q_{n} \in P_{n}} {‖ f - Q_{n} ‖}_{2}$ We also denote this by $P_{n} = \underset{Q_{n} \in P_{n}}{argmin} {‖ f - Q_{n} ‖}_{2}$ .

In order to calculate such a solution we need the following important result.

A necessary and sufficient condition for $P_{n} \in P_{n}$ being a best least square approximation of $f \in E$ is: $\forall Q_{n} \in P_{n}$ , $⟨ f - P_{n}, Q_{n} ⟩ = 0.$

This condition can be interpreted geometrically. Among all the polynoms $Q_{n}$ of $P_{n}$ the solution $P_{n}$ is obtained as the orthogonal projection of $f$ on the vector space $P_{n}$ which is included in $E$ .

$(\Rightarrow)$ Let $P_{n}$ be the best least square approximation of $f$ . We use a reductio ad absurdum argument (raisonnement par l'absurde). We assume there exists a polynom $Q \in P_{n}$ such that $⟨ f - P_{n}, Q ⟩ = α \neq 0$ . We set $R = P_{n} + λ Q$ . If $λ \in] 0, \frac{2 α}{{‖ Q ‖}_{2}^{2}} [$ then we have $\begin{array}{r} {‖ f - R ‖}_{2}^{2} = {‖ (f - P_{n}) - λ Q ‖}_{2}^{2} = {‖ (f - P_{n}) ‖}_{2}^{2} - 2 λ α + λ^{2} {‖ Q ‖}_{2}^{2} < {‖ f - P_{n} ‖}_{2}^{2} \end{array}$ By choosing $λ = \frac{α}{{‖ Q ‖}_{2}^{2}}$ we obtain a contradiction.
$(\Leftarrow)$ Let $P_{n} \in P_{n}$ such that for any $Q \in P_{n}$ , $⟨ f - P_{n}, Q_{n} ⟩ = 0$ . As $⟨ f - P_{n}, Q - P_{n} ⟩ = 0$ because $Q - P_{n} \in P_{n}$ , we then have ${‖ f - Q ‖}_{2}^{2} = ⟨ (f - P_{n}) - (Q - P_{n}), (f - P_{n}) - (Q - P_{n}) ⟩ = {‖ f - P_{n} ‖}_{2}^{2} + {‖ Q - P_{n} ‖}_{2}^{2}$ Consequently ${‖ f - P_{n} ‖}_{2} \leq {‖ f - Q ‖}_{2}$ and $P_{n}$ is the best least square approximation of $f$ .

Moreover, this projection is unique.

The best least square approximation $P_{n} \in P_{n}$ of $f \in E$ is unique.

Let $g_{1}$ and $g_{2}$ be two best least square approximation belonging to $P_{n}$ of $f$ . It yields $\begin{aligned} {‖ g_{1} - g_{2} ‖}_{2}^{2} & = ⟨ g_{1} - f + f - g_{2}, g_{1} - g_{2} ⟩ \\ = ⟨ g_{1} - f, g_{1} - g_{2} ⟩ + ⟨ f - g_{2}, g_{1} - g_{2} ⟩ \end{aligned}$ As $g_{1} - g_{2}$ belongs to $P_{n}$ we then have by the previous Theorem that the scalar products are equals to zero. Consequently $g_{1} = g_{2} .$

Exercice [PY-3]

We want to calculate $P_{2}$ , the best least square polynomial approximation of $f (x) = | x |$ on $[- 1, 1]$ . For this, we consider the following scalar product: $⟨ f, g ⟩ = \int_{a}^{b} f (x) g (x) d x$ .

By using the necessary and sufficient condition, show that the solution $P_{2} (x) = a_{0} + a_{1} x + a_{2} x^{2}$ can be obtained by solving the following linear system: $[\begin{array}{ccc} ⟨ 1, 1 ⟩ & ⟨ 1, x ⟩ & ⟨ 1, x^{2} ⟩ \\ ⟨ x, 1 ⟩ & ⟨ x, x ⟩ & ⟨ x, x^{2} ⟩ \\ ⟨ x^{2}, 1 ⟩ & ⟨ x^{2}, x ⟩ & ⟨ x^{2}, x^{2} ⟩ \end{array}] [\begin{matrix} a_{0} \\ a_{1} \\ a_{2} \end{matrix}] = [\begin{matrix} ⟨ f, 1 ⟩ \\ ⟨ f, x ⟩ \\ ⟨ f, x^{2} ⟩ \end{matrix}]$
Calculate all the coefficients of the previous system and solve it with python.
Draw the function and polynom $P_{2}$ .
If the scalar product becomes $⟨ f, g ⟩ = \int_{a}^{b} w (x) f (x) g (x) d x$ then the we have to solve the following system: $\underset{j = 0}{\sum^{p}} (\int_{- 1}^{1} w (x) x^{j + k} d x) a_{j} = \int_{- 1}^{1} w (x) f (x) x^{k} d x$ Calculate the new solution when $w (x) = \frac{1}{\sqrt{1 - x^{2}}}$ on $] - 1, 1 [$ .

            
import numpy as np
import time
from numpy import linalg as LA
import matplotlib.pyplot as plt


def f(x):
    return np.abs(x)

def Pol(x,coeff):
    s = 0.0
    # TO DO
    return s

# Define the matrix
def Mat(X):
    n=len(X)-1
    A=np.zeros((n+1,n+1))
    # TO DO
    return A

# degree of the polynom
n=2

#Plot function f
t=np.linspace(a,b,100)
p2=plt.plot(t,f(t),'k')

#TO DO: solve the linear system

# Plot the polynom
p3=plt.plot(t,Pol(t,coeff),'r')
        
    

                                            
import numpy as np
from numpy import linalg as LA
import matplotlib.pyplot as plt

# Interpolation polynomiale
print("Interpolation polynomiale")
def f(x):
    return np.abs(x)
    #return (x**2-4*x+4)*(x**2+4*x+2) 

# nombre de point d'interpolation
p=5
# degré du polynôme
n=p-1
# création de la liste des données à interpoler
a=-1.0
b= 1.0

# Horner Algorithm
def Pol(x,coeff):
    p=len(coeff)
    s=coeff[p-1]
    for i in range(p-1):
        s=x*s+coeff[p-2-i]
    return s


X=np.linspace(a,b,p)
Y=f(X)

t=np.linspace(a,b,1000)
plt.ylim(-1,2)
p1=plt.plot(X,Y,'bo')


# Création de la matrice d'interpolation
A=np.eye(n+1)
for i in range(n+1):
    for j in range(n+1):
        A[i,j]=X[i]**j    

# Resolution de AX=Y
coeff=LA.solve(A,Y)
print("coefficients = ",coeff)

p2=plt.plot(t,f(t),'k')
p3=plt.plot(t,Pol(t,coeff),'r')
plt.show()
                        
                    
                

References

C. Eckart, G. Young, The approximation of one matrix by another of lower rank. Psychometrika, Volume 1, 1936, Pages 211–8. doi:10.1007/BF02288367
C. de Boor (1978) A Practical Guide to Splines. Springer-Verlag.
G.D. Knott (2000) Interpolating Cubic Splines. Birkhäuser.
H.J. Nussbaumer (1981) Fast Fourier Transform and Convolution Algorithms. Springer-Verlag.
H. Späth (1995) One Dimensional Spline Interpolation. AK Peters.