Least Squares Problem

Example : \(C^{14}\) Datation

Radioactive carbon \(^{14}C\) is produced in the atmosphere by the effect of cosmic rays on atmospheric nitrogen. It is oxidized to \(^{14}CO_{2}\) and absorbed in this form by living organisms. So, living organisms contain a certain percentage of radioactive carbon relative to \(^{12}C\) and \(^{13}C\) which are stable. We suppose that carbon production \(^{14}C\) is constant over the last few millennia.

It is also assumed that, when an organism dies, its exchanges with the atmosphere cease, and that radioactivity due to carbon to carbon \(^{14}C\) decreases according to the following exponential law:

\[ A(t,A_0,\lambda)=A_{0}e^{-\lambda t}. \]

The analysis of the trunks (wood is dead tissue) of old trees Sequoia gigantea and {} furnishes us~:

its age \(t\) in year
its radioactivity \(A\)

\[ \begin{array}{||c|ccccccc||}\hline\hline t_i & 500 & 1000 & 2000 & 3000 & 4000 & 5000 & 6300 \\ \hline A_i & 14.5 & 13.5 & 12.0 & 10.8 & 9.9 & 8.9 & 8.0\\ \hline\hline \end{array} \]

We want to find the values of the parameters \(A_0\) and \(\lambda\), so that the function \(A{(t,A_0,\lambda)}\) is “near” the data :

\[(P) \left\{ \begin{array}{l} \displaystyle Min\; f(\beta) = \frac{1}{2}\|r(\beta\|^2 = f(A_0,\lambda)=\frac{1}{2}\sum_{i=1}^{n}(A_i-A_0e^{-\lambda t_i})^2\\ \beta=(A_0,\lambda) \in {\R}^2. \end{array} \right. \]

Solve this problem by using :

The Newton algorithm
The Gauss-Newton algorithm :

\[({P}_k)\left\{\begin{array}{l} Min\;\;f_k(s)=\frac{1}{2}\|r(\beta^{(k)})+J_r(\beta^{(k)})s\|^2\\ s\in \R^p, \end{array}\right.\] where \(s = \beta - \beta^{(k)}\) abd \(J_r(\beta)\) is the Jacobian matrix of \(r\) in \(\beta\).

Data, Residuals and \(f\) functions

Code

using LinearAlgebra
using Plots
Ti = [ 500, 1000, 2000, 3000, 4000, 5000, 6300];
Ai = [14.5, 13.5, 12.0, 10.8,  9.9,  8.9,  8.0];
n = length(Ti)
Data_C14 = [Ti Ai];
println(Data_C14)

# Initialisation
beta0 = [10; 0.0001];   # Newton, Gauus-Newton, fminunc et leastsq converge
# beta0 = [15; 0.001];    # Newton, Gauss-Newton, fminunc et leastsq divergent
# beta0 = [15; 0.0005];   # Newton diverge, Gauss-Newton, fminunc et leastsq convergent
# beta0 = [10; 0.0005];   # Gauss-Newton converge

# Initial model
#----------------------------------
xmin = 9; xmax = 20;
xx = range(xmin, stop=xmax, length=100);

ymin = -0.0001; ymax = 0.0005;
yy = range(ymin, stop=ymax, length=100);

# Residual function
function r(β,data)
    A₀ = β[1]
    λ = β[2]
    return data[:,2]-A₀*exp.(-λ*data[:,1])
end
# Plot of the function

println("r(beta0,Data_C14) = ", r(beta0,Data_C14))
X = [ones(n) Ti]
f(β) = 0.5*norm(r(β,Data_C14))^2
f_contour(A₀,λ) = f([A₀,λ])
z = @. f_contour(xx', yy)

p1 = plot()
contour!(p1,xx,yy,z,levels=100,cbar=false,color=:turbo)

p2 = scatter(Ti,Ai,title="Data C14")
T = range(0,stop=6500,length=100);
A = beta0[1]*exp.(-beta0[2]*T);
plot!(p2,T,A)

[500.0 14.5; 1000.0 13.5; 2000.0 12.0; 3000.0 10.8; 4000.0 9.9; 5000.0 8.9; 6300.0 8.0]
r(beta0,Data_C14) = [4.987705754992859, 4.451625819640405, 3.812692469220181, 3.391817793182822, 3.1967995396436066, 2.8346934028736666, 2.674081989931028]

With the Gauss-Newton Algorithm

Code

# solve by Gauss-Newton Algorithm
using Plots
include("assets/julia/MyOptims.jl")
nb_levels = 5
res(β) = r(β,Data_C14)
for nbit in 1:nb_levels
    βsol, flag, fsol, ∇f_xsol , nb_iter  = algo_Gauss_Newton(res,beta0,nbit_max=nbit)
 #   A = βsol[1]*exp.(-βsol[2]*T);
    plot!(p2,T,βsol[1]*exp.(-βsol[2]*T))
    scatter!(p1,[βsol[1]],[βsol[2]])
end

plot(p1,p2,legend=false)

With Newton Algorithm

Code

# solve by Newton Algorithm

nb_levels = 5
for nbit in 1:nb_levels
    βsol, flag, fsol, ∇f_xsol , nb_iter  = algo_Newton(f,beta0,nbit_max=nbit)
 #   A = βsol[1]*exp.(-βsol[2]*T);
    plot!(p2,T,βsol[1]*exp.(-βsol[2]*T))
    scatter!(p1,[βsol[1]],[βsol[2]])
end

plot(p1,p2,legend=false)