Explaining how the algorithm works and the concepts involved is quite difficult and so I shall do it slowly and section by section. I shall refer to the reader as Alice as in the famous character in the Novel "Alice through the Looking Glass". The reason for this is because I need to take the reader into the Looking Glass Universe where things are view differently

# Assumptions

In this article I assume that You know

- How to add a package to Julia
- What an XY graph is
- What a Matrix is
- The basic idea of what an inverse of a Matrix is
- The basic idea of what a transpose of a Matrix is
- Basic maths. Addition Subtration. Multiplication. Division. Power.

In this article, I will show you how to perform Least Square Fitting to obtain the parameters of a non-linear model on experimental data.

# Step 0 Talk abut the model

In this article we shall model the experimental data using the non-linear model

Y = A exp(B * X) + C

where X,Y are the measurement values from the experiment and

A,B,C are the unknown constants.

The unknown constants are also refer to as the parameters

We generalise this by saying

Y = f(x) and

f(x) = A exp(B * x) + C

Now this function **f** is **very important** and it represent the relationship between X and Y. That is to say: Given X, what is Y?

Later on in the article we will use the function **f** to get a "Looking Glass" function **g** by stepping through the Looking Glass. Naturally the function **g** is highly related to the function **f** . We shall talk about this later.

# Step 1 Obtaining experimental data

We could perform the experiment and obtain lots of values of X and Y but for this article we shall generate the values using a deterministic random number generator.

Here is the source code for generating the values of X and Y

```
using Distributions, StableRNGs
#=
We want to generate Float64 data deterministically
=#
function generatedata(formula; input=[Float64(_) for _ in 0:(8-1)], seed=1234, digits=0,
preamble="rawdata = ", noisedist=Normal(0.0, 0.0), verbose=true)
vprint(verbose, str) = print(verbose == true ? str : "")
local myrng = StableRNG(seed)
local array = formula.(input) + rand(myrng, noisedist, length(input))
if digits > 0
array = round.(array, digits=digits)
end
vprint(verbose, preamble)
# Now print the begining array symbol '['
vprint(verbose, "[")
local counter = 0
local firstitem = true
for element in array
if firstitem == true
firstitem = false
else
vprint(verbose, ",")
end
counter += 1
if counter % 5 == 1
vprint(verbose, "\n ")
end
vprint(verbose, element)
end
vprint(verbose, "]\n")
return array
end
A_actual = 1.5
B_actual = -0.25
C_actual = 3.5
f(x) = A_actual * exp(B_actual * x) + C_actual
numofsamples = 400+1
X = [Float64(_) for _ in 0:0.01:(numofsamples-1)/100]
Y = generatedata(f, input=X, digits=5, noisedist=Normal(0.0, 0.0005), preamble="Y = ",verbose=false)
println("\n\nX = ", X, "\nY = ", round.(Y,digits=5), "\n\n")
```

with the output

```
X = [0.0, 0.01, 0.02, 0.03, 0.04, 0.05, 0.06, 0.07, 0.08, 0.09, 0.1, 0.11, 0.12, 0.13, 0.14, 0.15, 0.16, 0.17, 0.18, 0.19, 0.2, 0.21, 0.22, 0.23, 0.24, 0.25, 0.26, 0.27, 0.28, 0.29, 0.3, 0.31, 0.32, 0.33, 0.34, 0.35, 0.36, 0.37, 0.38, 0.39, 0.4, 0.41, 0.42, 0.43, 0.44, 0.45, 0.46, 0.47, 0.48, 0.49, 0.5, 0.51, 0.52, 0.53, 0.54, 0.55, 0.56, 0.57, 0.58, 0.59, 0.6, 0.61, 0.62, 0.63, 0.64, 0.65, 0.66, 0.67, 0.68, 0.69, 0.7, 0.71, 0.72, 0.73, 0.74, 0.75, 0.76, 0.77, 0.78, 0.79, 0.8, 0.81, 0.82, 0.83, 0.84, 0.85, 0.86, 0.87, 0.88, 0.89, 0.9, 0.91, 0.92, 0.93, 0.94, 0.95, 0.96, 0.97, 0.98, 0.99, 1.0, 1.01, 1.02, 1.03, 1.04, 1.05, 1.06, 1.07, 1.08, 1.09, 1.1, 1.11, 1.12, 1.13, 1.14, 1.15, 1.16, 1.17, 1.18, 1.19, 1.2, 1.21, 1.22, 1.23, 1.24, 1.25, 1.26, 1.27, 1.28, 1.29, 1.3, 1.31, 1.32, 1.33, 1.34, 1.35, 1.36, 1.37, 1.38, 1.39, 1.4, 1.41, 1.42, 1.43, 1.44, 1.45, 1.46, 1.47, 1.48, 1.49, 1.5, 1.51, 1.52, 1.53, 1.54, 1.55, 1.56, 1.57, 1.58, 1.59, 1.6, 1.61, 1.62, 1.63, 1.64, 1.65, 1.66, 1.67, 1.68, 1.69, 1.7, 1.71, 1.72, 1.73, 1.74, 1.75, 1.76, 1.77, 1.78, 1.79, 1.8, 1.81, 1.82, 1.83, 1.84, 1.85, 1.86, 1.87, 1.88, 1.89, 1.9, 1.91, 1.92, 1.93, 1.94, 1.95, 1.96, 1.97, 1.98, 1.99, 2.0, 2.01, 2.02, 2.03, 2.04, 2.05, 2.06, 2.07, 2.08, 2.09, 2.1, 2.11, 2.12, 2.13, 2.14, 2.15, 2.16, 2.17, 2.18, 2.19, 2.2, 2.21, 2.22, 2.23, 2.24, 2.25, 2.26, 2.27, 2.28, 2.29, 2.3, 2.31, 2.32, 2.33, 2.34, 2.35, 2.36, 2.37, 2.38, 2.39, 2.4, 2.41, 2.42, 2.43, 2.44, 2.45, 2.46, 2.47, 2.48, 2.49, 2.5, 2.51, 2.52, 2.53, 2.54, 2.55, 2.56, 2.57, 2.58, 2.59, 2.6, 2.61, 2.62, 2.63, 2.64, 2.65, 2.66, 2.67, 2.68, 2.69, 2.7, 2.71, 2.72, 2.73, 2.74, 2.75, 2.76, 2.77, 2.78, 2.79, 2.8, 2.81, 2.82, 2.83, 2.84, 2.85, 2.86, 2.87, 2.88, 2.89, 2.9, 2.91, 2.92, 2.93, 2.94, 2.95, 2.96, 2.97, 2.98, 2.99, 3.0, 3.01, 3.02, 3.03, 3.04, 3.05, 3.06, 3.07, 3.08, 3.09, 3.1, 3.11, 3.12, 3.13, 3.14, 3.15, 3.16, 3.17, 3.18, 3.19, 3.2, 3.21, 3.22, 3.23, 3.24, 3.25, 3.26, 3.27, 3.28, 3.29, 3.3, 3.31, 3.32, 3.33, 3.34, 3.35, 3.36, 3.37, 3.38, 3.39, 3.4, 3.41, 3.42, 3.43, 3.44, 3.45, 3.46, 3.47, 3.48, 3.49, 3.5, 3.51, 3.52, 3.53, 3.54, 3.55, 3.56, 3.57, 3.58, 3.59, 3.6, 3.61, 3.62, 3.63, 3.64, 3.65, 3.66, 3.67, 3.68, 3.69, 3.7, 3.71, 3.72, 3.73, 3.74, 3.75, 3.76, 3.77, 3.78, 3.79, 3.8, 3.81, 3.82, 3.83, 3.84, 3.85, 3.86, 3.87, 3.88, 3.89, 3.9, 3.91, 3.92, 3.93, 3.94, 3.95, 3.96, 3.97, 3.98, 3.99, 4.0]
Y = [5.00026, 4.99671, 4.99168, 4.98814, 4.98473, 4.98102, 4.97771, 4.97393, 4.96957, 4.96567, 4.96252, 4.96015, 4.9556, 4.95183, 4.94868, 4.94439, 4.94139, 4.93792, 4.93397, 4.93044, 4.92662, 4.92324, 4.91925, 4.91613, 4.91287, 4.90931, 4.90609, 4.90147, 4.89842, 4.89508, 4.89157, 4.8876, 4.88476, 4.88175, 4.87781, 4.87462, 4.87074, 4.86746, 4.86411, 4.86027, 4.8575, 4.85359, 4.84989, 4.84682, 4.84404, 4.83992, 4.83673, 4.83395, 4.83077, 4.82657, 4.82314, 4.82066, 4.81735, 4.81311, 4.80956, 4.80774, 4.80396, 4.80025, 4.79692, 4.79453, 4.79056, 4.78834, 4.78521, 4.78146, 4.77902, 4.77504, 4.77291, 4.76972, 4.7647, 4.76251, 4.75913, 4.75643, 4.75262, 4.7496, 4.74612, 4.74393, 4.74016, 4.73748, 4.73441, 4.73105, 4.72772, 4.72505, 4.72198, 4.71788, 4.71669, 4.71308, 4.71046, 4.70721, 4.70349, 4.70032, 4.69721, 4.69474, 4.69145, 4.68943, 4.68538, 4.68237, 4.68024, 4.67703, 4.67413, 4.67229, 4.66888, 4.66589, 4.66264, 4.65933, 4.65577, 4.65449, 4.65059, 4.64823, 4.6461, 4.6431, 4.6392, 4.63667, 4.63245, 4.63044, 4.62859, 4.62477, 4.62197, 4.62009, 4.61724, 4.61425, 4.6117, 4.6092, 4.60616, 4.60314, 4.59975, 4.5978, 4.59554, 4.59203, 4.58831, 4.58639, 4.58263, 4.58192, 4.57885, 4.5761, 4.57306, 4.57052, 4.56772, 4.56495, 4.56203, 4.55984, 4.55683, 4.55441, 4.55088, 4.54908, 4.54606, 4.54433, 4.54214, 4.53836, 4.5362, 4.53241, 4.53141, 4.52867, 4.52551, 4.52338, 4.52062, 4.51883, 4.51532, 4.51249, 4.51068, 4.50698, 4.50629, 4.50336, 4.50098, 4.49811, 4.49508, 4.49317, 4.48983, 4.48791, 4.48504, 4.48337, 4.48136, 4.4778, 4.47559, 4.47316, 4.4703, 4.46867, 4.46528, 4.46452, 4.45991, 4.45953, 4.45685, 4.45328, 4.45157, 4.44927, 4.44621, 4.44493, 4.4421, 4.43922, 4.43794, 4.43555, 4.43322, 4.43068, 4.42809, 4.4264, 4.42316, 4.42128, 4.41914, 4.41739, 4.41406, 4.41157, 4.40969, 4.4076, 4.40453, 4.40268, 4.40096, 4.39855, 4.39654, 4.39399, 4.39193, 4.3892, 4.38635, 4.38445, 4.38312, 4.38146, 4.37867, 4.37581, 4.37449, 4.37264, 4.3692, 4.36779, 4.36556, 4.36366, 4.36077, 4.3589, 4.3578, 4.35527, 4.35167, 4.34997, 4.3491, 4.3457, 4.34382, 4.34193, 4.3396, 4.33865, 4.33591, 4.33407, 4.33229, 4.32937, 4.32821, 4.3259, 4.32237, 4.32102, 4.31951, 4.31657, 4.31483, 4.31291, 4.31125, 4.30921, 4.30641, 4.30484, 4.3028, 4.301, 4.29897, 4.29732, 4.29448, 4.29281, 4.28975, 4.28915, 4.2868, 4.2858, 4.28295, 4.28058, 4.28009, 4.27766, 4.27502, 4.27358, 4.27174, 4.2687, 4.26786, 4.26582, 4.26392, 4.26183, 4.25931, 4.25732, 4.25707, 4.25473, 4.25133, 4.24999, 4.24787, 4.24655, 4.24507, 4.24296, 4.24149, 4.23917, 4.23702, 4.2357, 4.23442, 4.23145, 4.23031, 4.22895, 4.22615, 4.22453, 4.22319, 4.22086, 4.21921, 4.21783, 4.21609, 4.21384, 4.21163, 4.21076, 4.20879, 4.2065, 4.20493, 4.20258, 4.20095, 4.19956, 4.19858, 4.19528, 4.19366, 4.19288, 4.1901, 4.18912, 4.1877, 4.18586, 4.18422, 4.18253, 4.18112, 4.17955, 4.17818, 4.17609, 4.17396, 4.17281, 4.17105, 4.16854, 4.16699, 4.16667, 4.16367, 4.16227, 4.16073, 4.15925, 4.15801, 4.15625, 4.15331, 4.15193, 4.15034, 4.14909, 4.14783, 4.14581, 4.14392, 4.14243, 4.1407, 4.13978, 4.13744, 4.13692, 4.13441, 4.13327, 4.13116, 4.13048, 4.12872, 4.12701, 4.1243, 4.12402, 4.12258, 4.11982, 4.11875, 4.11768, 4.11669, 4.1146, 4.11309, 4.11081, 4.10998, 4.10841, 4.10661, 4.10511, 4.10405, 4.10142, 4.10056, 4.09889, 4.09851, 4.09635, 4.09528, 4.09188, 4.09195, 4.0895, 4.08858, 4.08764, 4.0859, 4.08428, 4.08316, 4.08134, 4.07975, 4.07809, 4.07754, 4.07452, 4.0756, 4.07306, 4.07144, 4.06975, 4.06911, 4.06743, 4.06596, 4.06462, 4.06311, 4.06288, 4.05962, 4.05898, 4.05684, 4.05688, 4.05526, 4.05294, 4.05135]
```

# Step 2 Doing the Least Square Fitting on the experimental data

We show you the source code for the algorithmn and the output first and later in the article we shall explain how the algorithmn works.

```
"""
Provide the solution to the nonlinear Least Square
"""
function nonlinearLSQ(xarray, yarray; maxiter=8, startpos=[ 1.0 , -0.1, 1.0 ], stepsize=1.0, verbose=true)
len_xarray = length(xarray)
len_yarray = length(yarray)
if len_xarray != len_yarray
println("Error in function nonlinearLSQ: length of xarray != length of yarray")
return nothing
end
# Parameter array P = [ A , B , C ]
P = startpos # Initial Starting Position
for iteration = 1:maxiter
R = zeros(len_yarray) # a len_yarray size array
# Let's create the J Matrix (Jacobian Matrix)
J = zeros(len_yarray, 3) # a len_yarray*3 matrix
for i = 1:len_yarray
J[i, 1] = exp(P[2] * X[i])
J[i, 2] = P[1] * exp(P[2] * X[i]) * X[i]
J[i, 3] = 1.0
R[i] = yarray[i] - ( P[1] * exp(P[2] * xarray[i]) + P[3] )
end
# Dont forget to reshape values from
# a 1D vector R into a len*1 column matrix RR
local RR = reshape(R, len_yarray, 1) # a len_yarray * 1 Matrix
DeltaP = inv(J' * J) * J' * RR
# At this point the varible DeltaP is a len_yarray * 1 Matrix
# We turn it into an 1 dimensional vector using vec() function
P = P + stepsize .* vec( DeltaP )
if verbose == true
println("Iteration is ",iteration)
println(" DeltaP is ",round.(DeltaP,digits=5))
println(" P is ",round.(P,digits=5))
end
end
println()
return P
end
solution = nonlinearLSQ(X, Y, maxiter=8, startpos=[ 1.0 , -0.1, 1.0 ])
println("Trying with startpos [ 1.0 , -0.1 , 1.0 ] stepsize=1.0 maxiter=8")
println("The estimate for A is ", round.(solution[1],digits=5))
println("The estimate for B is ", round.(solution[2],digits=5))
println("The estimate for C is ", round.(solution[3],digits=5))
```

With the output

```
Iteration is 1
DeltaP is [-1.58274; -0.42322; 4.57972;;]
P is [-0.58274, -0.52322, 5.57972]
Iteration is 2
DeltaP is [1.82102; -0.53788; -1.82729;;]
P is [1.23828, -1.0611, 3.75242]
Iteration is 3
DeltaP is [-0.34721; 0.78483; 0.28928;;]
P is [0.89107, -0.27627, 4.0417]
Iteration is 4
DeltaP is [0.601; 0.04232; -0.53397;;]
P is [1.49208, -0.23395, 3.50773]
Iteration is 5
DeltaP is [0.00431; -0.01652; -0.00425;;]
P is [1.49638, -0.25047, 3.50349]
Iteration is 6
DeltaP is [0.00429; 0.00069; -0.00425;;]
P is [1.50067, -0.24979, 3.49923]
Iteration is 7
DeltaP is [1.0e-5; -0.0; -1.0e-5;;]
P is [1.50068, -0.24979, 3.49923]
Iteration is 8
DeltaP is [0.0; 0.0; -0.0;;]
P is [1.50068, -0.24979, 3.49923]
Trying with startpos [ 1.0 , -0.1 , 1.0 ] stepsize=1.0 maxiter=8
The estimate for A is 1.50068
The estimate for B is -0.24979
The estimate for C is 3.49923
```

# Step 3 Discussing how Jacobian Matrix is constructed in the code.

In the source code, it was mentioned that you need to create the Jacobian Matrix. Now you may want to read the wikipedia article below on how to create the Jacobian Matrix but I found the wikipedia article to be confusing for newbies.

https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant

The reason the wikipedia article is confusing is that it talks about a function f but that f is referred to a particular generic function related to the issue at hand.

I have modified the content of the wikipedia article to **better** reflect our problem and it looks like this.

Now time for some explaination. The Matrix **p** is a 3x1 parameter matrix. It consists of

Now the function **g** on the otherhand is a family of functions.

Because we have 401 data points from the experiment. Namely

Now Alice, we come back to the function **f** which we mentioned at the start of the article. As you recall, we had the non-linear model

Y = A exp(B * X) + C

where X,Y are the measurement values from the experiment and

A,B,C are the unknown constants.

Which we generalise by saying

Y = f(x) and

f(x) = A exp(B * x) + C

Now Alice, understand that this is **TRUE** before we perform the experiment. That is to say X and Y are variables while A, B and C are the unknown constants.

That is to say

X can be any value(s)

Y can be any value(s)

but

A is a fixed value (but currently unknown)

B is a fixed value (but currently unknown)

C is a fixed value (but currently unknown)

And the **WHOLE POINT** of the experiment is to determine the values of A, B and C

Now it is time to enter the Looking Glass by performing the experiment. Once we had perform the experiment, we are in the **MAGICAL REALM** of the **LOOKING GLASS** and things change:

X and Y are now a list of **FIXED** value pairs

x1 , y1

x2 , y2

**.**

**.**

**.**

x401 , y401

We represent this by saying **x** and **y** is a 401x1 column matrix

And we represent it by

**y** = **g**(**p**)

Where **g** is a family of functions consists of

**g** is a 1x401 matrix of

g1(A,B,C) = A exp(B * x1) + C

g2(A,B,C) = A exp(B * x2) + C

g3(A,B,C) = A exp(B * x3) + C

**.**

**.**

**.**

g401(A,B,C) = A exp(B * x401) + C

Which in this case, we can say g(A,B,C) is a function related to the function of f(x)

Now when we put in the values of x1 , x2 , x3 , ... , x401 from the experimental results, we get

g1(A,B,C) = A exp(B * 0.0) + C

g2(A,B,C) = A exp(B * 0.01) + C

g3(A,B,C) = A exp(B * 0.02) + C

**.**

**.**

**.**

g401(A,B,C) = A exp(B * 4.0) + C

It is important to know WHAT HAS CHANGED

That is to say X and Y are constants while A, B and C are now variables

That is to say

X is FIXED (by the values obtained from the experiment)

Y is FIXED (by the values obtained from the experiment)

but

A can be any value

B can be any value

C can be any value

We need to do this because we wanted to differentiate the variables A , B and C and we cannot (technically) differentiate a constant.

So back to generating the Jacobian Matrix

We need to find the partial derivative of g1 in respect to A

in particular

g1(A,B,C) = A exp(B * x1) + C

There are 3 ways

- Using Mathematica
- Using Wolfram Alpha
- Doing it by hand

## Using Mathematica

## Using Wolfram Alpha

using url https://www.wolframalpha.com

## Doing it by hand

Too lazy to do it by hand, Here are the basic steps

The equivalent part of the source code is

```
J[i, 1] = exp(P[2] * X[i])
```

Remember P[2] is the variable B

## Partial derivative with respect to B

Next We need to find the partial derivative of g1 in respect to B

in particular

g1(A,B,C) = A exp(B * x1) + C

The answer is A * exp(B * x1) * x1

Which is reflected in the source code

```
J[i, 2] = P[1] * exp(P[2] * X[i]) * X[i]
```

## Partial derivative with respect to C

Finally We need to find the partial derivative of g1 in respect to C

in particular

g1(A,B,C) = A exp(B * x1) + C

The answer is 1

Which is reflected in the source code

```
J[i, 3] = 1.0
```

# Step 4 Discussing how the algorithm works. Multiple explainations. Pick one that you like

## How does the algorithmn work?

There are many explanation of the algorithm works. Here are some of the explainations

# Step 5 Discussing the effect of the step size

Our starting pos is

solution = nonlinearLSQ(X, Y, maxiter=8, startpos=[ 1.0 , -0.1, 1.0 ])

So why do we start with startpos = [ 1.0 , -0.1 , 1.0 ]? Because I knew the answer is 1.5 , -0.25 , 3.5 . I know the answer is that because I builted up the raw data.

But what if we do not know what the answer is?

Then I woud have started the program with the startpos of [ 1.0 , -1.0 , 1.0 ] being a generic starting point. I tend to avoid values of 0.0

So what happens when we start with that startpos instead?

## Top comments (0)