Computational Methods
for Physicists and Materials Engineers
5
Least squares fitting
and machine learning (brief introduction)
Recall Lecture 3 : Linear regression
Problem: We want to find a function y = y(x(1), ∙∙∙, x(N)) E.g. y = lifetime of a material before failure
x(1) = service temperature
x(2) = humidity of the environment
x(3), x(4) = hydrostatic and deviatoric stresses …
The simplest guess is the linear relationship
y w(0) w(1)x(1) w(2)x(2) w(N)x(N)
(ǂ)
Now, the problem is to determine {w(i)} Collect a set of data by measurement:
yy whenx(1) x(1),x(2) x(2),,x(N) x(N) 1111
yy whenx(1) x(1),x(2) x(2),,x(N) x(N) 2222
yy whenx(1) x(1),x(2) x(2),,x(N) x(N)
MMMM Determine {w(i)} s.t. all the data are best described by Eq. (ǂ)
Recall Lecture 3 : Linear regression
y 1x(1)x(N) w(0) 11 1
yy2 ?1 x(1) x(N) w(1) 22
Xw
y 1 x(1) x(N) w(N)
MM1 M M M(N1) (N1)1 Least squares fitting: minimize the cost function J w.r.t. w
J(w)1 Xwy 2 1(Xwy)(Xwy)1w(XTX)wwXTy1yy 22222
Optimization problem: find w s.t. J(w) is minimized
J(w) 0 (XT X)w XT y 0 (XT X)w XT y
w
XT X w XT y Normal
(N1)(N1) (N1)1 (N1)1 equations Axb
Linear regression via QR decomposition (direct method)
Normal equations
XT Xw XT y
R
Recall QR decomposition
A=Q QTQ I
T 1 T XQR T T w(X X) X y(QR) (QR)1(QR) y
RT QT QR1 RT QT y RT R1 RT QT y
R1RTRTQT y R1QT y Rw QT y
Recall that R is an upper (right) triangular matrix. We can solve this system of equations simply be backward substitution
Linear regression via QR decomposition (direct method) XQRQ1R1
X = Q1 Q2
= Q1
R1
w= Q1T
y
R wQTy 11
0
R1
R1
Linear regression via QR decomposition (direct method)
# in file “regression.py”
import numpy as np
import scipy.linalg
def linear_regression_qr(X, y):
Implementation by Python
# add a column [1, …, 1]^T as 1st column of X
X = np.concatenate((np.ones((X.shape[0], 1)), X), axis=1)
# QR decomposition
Q, R = scipy.linalg.qr(X)
col_R = R.shape[1]
Q = Q[:, :col_R]; R = R[:col_R, :]
# backward substitution
b = np.matmul(np.transpose(Q), y)
N = b.shape[0]
w = np.empty(shape=b.shape)
w[N-1] = b[N-1] / R[N-1, N-1]
for m in range(N-2, -1, -1):
return w
w[m] = b[m]
for k in range(m+1, N):
w[m] -= R[m,k] * w[k]
w[m] = w[m] / R[m,m]
Linear regression via QR decomposition (direct method)
We can use the code to fit data by a plane in 3D
# in file “plane_fitting.py”
import numpy as np
from regression import linear_regression_qr
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
# generate data
x1 = 10 * np.random.rand(1000) – 5
x2 = 10 * np.random.rand(1000) – 5
y = 3 + x1 + 2 * x2 + np.random.normal(0, 1, x1.size)
# fitting
x1_col = np.reshape(x1, (x1.size, 1))
x2_col = np.reshape(x2, (x2.size, 1))
y_col = np.reshape(y, (y.size, 1))
X = np.concatenate((x1_col, x2_col), axis=1) w = linear_regression_qr(X, y_col)
x1_f = np.linspace(np.min(x1), np.max(x1), 1000, endpoint=True) x2_f = np.linspace(np.min(x2), np.max(x2), 1000, endpoint=True) X1_f, X2_f = np.meshgrid(x1_f, x2_f)
Y_f = w[0] + w[1] * X1_f + w[2] * X2_f
Linear regression via QR decomposition (direct method)
Plot original data and the plane obtained by fitting
# in file “plane_fitting.py”
fig, axs = plt.subplots(subplot_kw={‘projection’:’3d’}) axs.plot(x1, x2, y, ‘o’, markersize=2, color=’k’) axs.plot_surface(X1_f, X2_f, Y_f, linewidth=0, cmap=’rainbow’, alpha=0.5)
plt.show()
Linear regression via gradient descent (iterative method)
Least squares fitting: minimize the cost function J w.r.t. w J(w)1 Xwy2
2M
Optimization problem: find w s.t. J(w) is minimized
Guess the initial w0. Update w by iteration:
w(0) w(0) α J
w wαJ w(1)w(1)αJ
n1 n
w n1 n w(1) w w n
Gradient
n1 n w(0)
wwn w w n
w(N) w(N) α J n1 n w(N)
ww
n
Linear regression via gradient descent (iterative method) wn1 wn α J
How to understand gradient descent? Position of a ball: x and potential
energy: E(x). Force on the ball:
f E
x
If this ball moves along the force with velocity proportional to force,
finally the ball will stop at the bottom, i.e., the minimum of E E
xn1 xn αf(xn)xn αE
x xxn
w wwn Gradient
x
Linear regression via gradient descent (iterative method)
J(w)
Xw y
x(0)w(0) x(N)w(N) y 2
121M
2M
J 1
2
iii
w(j) M
M
iiii
i1
J α
M n1 n w(0) n Mi1
x(0)w(0) x(N)w(N) y x(0) ininii
w(0) w(0) α w(0) wwn
M J α
x(0)w(0) x(N)w(N) y x(1) ininii
w(1) w(1) α w(1)
n1 n w(1) n M w w n
i1 M
J α w(N) w(N) α w(N)
x(0)w(0) x(N)w(N) y x(N) ininii
n 1 n w ( N ) n M i 1
2M i1
x(0)w(0) x(N)w(N) y x( j)
ww
Vectorized form (for NumPy element‐wise operation)
n
w w αJ w αXTXw y n1 n w n M n
wwn
Linear regression via gradient descent (iterative method)
# in file “regression.py”
import numpy as np
import scipy.linalg
Implementation by Python
def linear_regression_gd(X, y, alpha=0.01, tol=1e-5, maxiter=10000000):
# add a column [1, …, 1]^T as 1st column of X
X = np.concatenate((np.ones((X.shape[0], 1)), X), axis=1) M, N = X.shape
w = np.zeros((N, 1)) # initial guess
is_converge = False
for i in range(maxiter + 1):
gradient = np.matmul(np.transpose(X), np.matmul(X, w) – y) if np.linalg.norm(gradient) < tol:
is_converge = True; break
w = w - (alpha / M) * gradient
if not is_converge: print("GD has not converged!") return w
# in file “plane_fitting.py”
# change the line “w = linear_regression_qr(X, y_col)” to w = linear_regression_gd(X, y_col)
Polynomial regression
Fit data to the following relationship (curve fitting): y(x)w(0) w(1)xw(2)x2 w(N)xN
There is one variable or feature, x. But this is a nonlinear function of x Convert the single‐variable nonlinear function to a multi‐variable linear
function. Let
x(1) x, x(2) x2, , x(N) xN
y w(0) w(1)x(1) w(2)x(2) w(N)x(N)
Then, we can simply apply linear regression
We can also fit data to the relationship (surface fitting):
y(x(1) ,x(2))w(0) w(1)x(1) w(2)x(2) w(3)(x(1))2 w(4)x(1)x(2) w(5)(x(2))2
Convert the multi‐variable nonlinear function to a multi‐variable linear
function. Let
x(3) (x(1) )2 , x(4) x(1) x(2) , x(5) (x(2) )2
Hexagonal close‐packed (hcp) crystal Lattice parameters: a, c
How can we calculate the equilibrium lattice parameters a0, c0?
Calculate the energy per atom under different (a, c) by first‐principles method: {Ei(ai, ci)}
Data:
c
Polynomial regression
Implementation by Python
x(1) x(2) E1 1a1 c1
E 1ac y 2 , X 2 2
E 1a c
a
y(x(1) ,x(2))w(0) w(1)x(1) w(2)x(2) w(3)(x(1))2 w(4)x(1)x(2) w(5)(x(2))2
MMM Fit the data to
# in file “regression.py”
import numpy as np
import scipy.linalg
Polynomial regression
Implementation by Python
def regression_surface(x1, x2, y, x1plot=None, x2plot=None):
# convert lists to numpy.arrays
x1 = np.array(x1); x2 = np.array(x2); y = np.array(y)
# convert to column vectors
M = y.size
x1 = np.reshape(x1, (M, 1))
x2 = np.reshape(x2, (M, 1))
y = np.reshape(y, (M, 1))
X = np.concatenate((x1, x2, x1**2, x1*x2, x2**2), axis=1) w = linear_regression_qr(X, y)
if x1plot is None or x2plot is None:
return w
else:
yplot = w[0] + w[1] * x1plot + w[2] * x2plot \
+ w[3] * x1plot**2 + w[4] * x1plot * x2plot \ + w[5] * x2plot**2
return yplot
# in file “surface_fitting.py”
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D from regression import regression_surface
def split_string_to_data(string): ...
lattparam_a = [];
lattparam_c = [];
energy = [];
file = open("hcp_fit.csv", 'r')
for line_index, line_str in enumerate(file):
if line_index == 0:
continue # skip the 1st line
Polynomial regression
Implementation by Python
# read data from “hcp_fit.csv”; see code in lecture 1
line_list = split_string_to_data(line_str) lattparam_a.append(float(line_list[0])) lattparam_c.append(float(line_list[1])) energy.append(float(line_list[2]))
file.close()
(continue last page)
Polynomial regression
Implementation by Python
a = np.linspace(np.min(lattparam_a), np.max(lattparam_a), 2000, endpoint=True)
c = np.linspace(np.min(lattparam_c), np.max(lattparam_c), 2000, endpoint=True)
A, C = np.meshgrid(a, c)
E = regression_surface(lattparam_a, lattparam_c, energy, A, C)
print(rf"Estimate of a is {A[np.where(E == np.min(E))]} Angstrom.")
print(rf"Estimate of c is {C[np.where(E == np.min(E))]} Angstrom.")
# plot the raw data and the surface obtained by fitting
fig, axs = plt.subplots(subplot_kw={'projection':'3d'})
axs.plot(lattparam_a, lattparam_c, energy, 'o', markersize=4, color='k')
axs.plot_surface(A, C, E, linewidth=0, cmap='rainbow', alpha=0.5) plt.tight_layout()
plt.show()
Polynomial regression
Implementation by Python
Fit data to the following relationship:
y f(x(1),,x(N);w)
Nonlinear regression
where f is a nonlinear equation, e.g. f(x; a, b) = tanh(ax + b) Again, least squares fitting is to minimize the cost function:
J(w)1 f(X;w)y21(fTf2yTfyTy) 2M 22M
1 M
2M (fi fi 2yi fi yiyi)
i1
Normal equation does not exist. So, we have to use gradient descent
J w(k)
M
f J 1 f T
(f y ) i (fy)
1
ii w(k) w M w
M
i1
wn1 wn α J
w wwn
import numpy as np
import scipy.stats
import matplotlib.pyplot as plt
Least squares fitting by Python
Linear regression, single‐variate (line fitting)
scipy.stats.linregress()
# generate raw data (in practice data are obtained by experiments)
num_points = 500
xmin = 0; xmax = 10
x = (xmax - xmin) * np.random.rand(num_points) + xmin y = 2 * x + 10
mean = 0; stddev = 2
y += np.random.normal(mean, stddev, num_points)
# plot raw data
fig, ax = plt.subplots()
ax.plot(x, y, 'x', color='#0000FF', label='Data') ax.set_xlim([xmin - 1, xmax + 1]) ax.set_ylim([5, 35])
fig.set_size_inches(5, 4)
plt.tight_layout()
plt.show()
Least squares fitting by Python
Linear regression, single‐variate (line fitting)
scipy.stats.linregress()
Least squares fitting by Python
Linear regression, single‐variate (line fitting)
scipy.stats.linregress()
res = scipy.stats.linregress(x, y)
print(f"Fitted slope: {res.slope}; error: {res.stderr}") print(f"Fitted intercept: {res.intercept}; error: {res.intercept_stderr}")
fig, ax = plt.subplots()
ax.plot(x, y, 'x', color='#0000FF', label='Data') ax.plot(x, res.intercept + res.slope * x, '-', linewidth=4, color='#A81551', label='Fitting')
ax.set_xlim([xmin-1, xmax+1])
ax.set_ylim([5, 35])
ax.set_xticks([0, 2, 4, 6, 8, 10])
ax.set_yticks([5, 10, 15, 20, 25, 30, 35]) ax.tick_params(axis='x', labelsize=10) ax.tick_params(axis='y', labelsize=10)
ax.grid(True)
ax.legend(fontsize=10, fancybox=False, framealpha=1.0) fig.set_size_inches(5,4)
plt.tight_layout()
plt.show()
Least squares fitting by Python
Linear regression, single‐variate (line fitting)
scipy.stats.linregress()
Least squares fitting by Python
Linear/Nonlinear regression, single‐variate (curve fitting)
Rational function
f(x)w(0) w(1)xw(2)x2 w(3)x3 1w(4)x w(5)x2 w(6)x3
import numpy as np
import scipy.optimize
import matplotlib.pyplot as plt
Least squares fitting by Python
Linear/Nonlinear regression, single‐variate (curve fitting)
scipy.optimize.curve_fit()
from file_handling import split_string_to_data
# read data from file
x_dat = []; y_dat = []
file = open("thermal_expansion.txt", 'r') for line_index, line_str in enumerate(file):
if line_index == 0:
continue # skip the 1st line
line_list = split_string_to_data(line_str)
# see lecture note 2 for function “split_string_to_data()” y_dat.append(float(line_list[0])) x_dat.append(float(line_list[1]))
file.close()
x_dat = np.array(x_dat); y_dat = np.array(y_dat) sort_ind = np.argsort(x_dat)
x_dat = x_dat[sort_ind]; y_dat = y_dat[sort_ind]
# define the fitting function
Least squares fitting by Python
Linear/Nonlinear regression, single‐variate (curve fitting)
scipy.optimize.curve_fit()
f(x)w(0) w(1)xw(2)x2 w(3)x3 1w(4)x w(5)x2 w(6)x3
def therm_expansion_func(x, w0, w1, w2, w3, w4, w5, w6): return (w0 + w1*x + w2*x**2 + w3*x**3) \
/ (1. + w4*x + w5*x**2 + w6*x**3)
popt, pcov = scipy.optimize.curve_fit(therm_expansion_func, x_dat, y_dat)
# x_dat and y_dat are arrays
# popt: optimized fitting parameters
# pcov: covariance matrix. diagonal element is variance of each fitting parameter
print(f"Fitting parameters are: \n{popt}")
print(f"Errors are: \n{np.sqrt(np.diag(pcov))}")
# error is standard deviation, i.e., sqrt of variance
# data points on fitting curve
fig, ax = plt.subplots()
Least squares fitting by Python
Linear/Nonlinear regression, single‐variate (curve fitting)
scipy.optimize.curve_fit()
x_plot = np.linspace(np.min(x_dat), np.max(x_dat), 1000, endpoint=True)
y_plot = therm_expansion_func(x_plot, *popt) # *popt: unpack popt
ax.plot(x_dat, y_dat, 'x', color='#0000FF', label='Data')
# plot raw data
ax.plot(x_plot, y_plot, '-', linewidth=2, color='#A81551', label='Fitting') # plot fitting curve
ax.set_xlim([0, 900])
ax.set_ylim([0, 22])
ax.set_xlabel('Temperature (K)', fontsize=12) ax.set_ylabel('Coefficient of thermal expansion', fontsize=12) ax.set_title('Thermal expansion of copper', fontsize=12) ax.grid(True)ax.legend(fontsize=10, fancybox=False, framealpha=1.0)
fig.set_size_inches(5,4)
plt.tight_layout()
plt.show()
Least squares fitting by Python
Linear/Nonlinear regression, single‐variate (curve fitting)
scipy.optimize.curve_fit()
Energy as a function of lattice parameters a and c for hcp crystals obtained by DFT calculation
Least squares fitting by Python
Linear/Nonlinear regression, multi‐variate (surface fitting)
scipy.optimize.curve_fit()
Quadratic function
E(a,c)w(0) w(1)aw(2)cw(3)a2 w(4)acw(5)c2
import numpy as np
import scipy.optimize
Least squares fitting by Python
Linear/Nonlinear regression, multi‐variate (surface fitting)
scipy.optimize.curve_fit()
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from file_handling import split_string_to_data
# read data from file
a_dat = []; c_dat = []; e_dat = []
file = open("hcp_fit.csv", 'r')
for line_index, line_str in enumerate(file):
if line_index == 0:
continue # skip the 1st line
line_list = split_string_to_data(line_str)
# see lecture note 2 for function “split_string_to_data()” a_dat.append(float(line_list[0])) c_dat.append(float(line_list[1])) e_dat.append(float(line_list[2]))
file.close()
Least squares fitting by Python
Linear/Nonlinear regression, multi‐variate (surface fitting)
scipy.optimize.curve_fit()
E(a,c)w(0) w(1)aw(2)cw(3)a2 w(4)acw(5)c2
# define the fitting function
def energy_lattparam(X, w0, w1, w2, w3, w4, w5): return w0 + w1*X[0,:] + w2*X[1,:] \
+ w3*X[0,:]**2 + w4*X[0,:]*X[1,:] + w5*X[1,:]**2
Xa_dat = np.array(a_dat).reshape((1, len(a_dat))) Xc_dat = np.array(c_dat).reshape((1, len(c_dat))) X_dat = np.concatenate((Xa_dat, Xc_dat), axis=0) # X_dat = [[ ... Xa_dat ... ],
# [ ... Xc_dat ... ]]
popt, pcov = scipy.optimize.curve_fit(energy_lattparam, X_dat, e_dat)
# Dimension of X_dat must be M x N. M is the number of variables (features). N is the number of input data (measurements) print(f"Fitting parameters are: \n{popt}")
print(f"Errors are: \n{np.sqrt(np.diag(pcov))}")
# data points on fitting curve
Least squares fitting by Python
Linear/Nonlinear regression, multi‐variate (surface fitting)
scipy.optimize.curve_fit()
a = np.linspace(np.min(a_dat), np.max(a_dat), 2000, endpoint=True) c = np.linspace(np.min(c_dat), np.max(c_dat), 2000, endpoint=True) A_mat, C_mat = np.meshgrid(a, c)
A_row = A_mat.reshape((1, -1)) # list matrix elements in a column C_row = C_mat.reshape((1, -1))
X_plot = np.concatenate((A_row, C_row), axis=0)
E_row = energy_lattparam(X_plot, *popt).reshape((1, -1))
# convert 3 row vectors: A_row, C_row, E_row to 2D meshgrid A_col = A_row.T; C_col = C_row.T; E_col = E_row.T
cols = np.unique(A_col).shape[0]
A = A_col.reshape(-1, cols)
C = C_col.reshape(-1, cols)
E = E_col.reshape(-1, cols)
fig, axs = plt.subplots(subplot_kw={'projection':'3d'})
axs.plot(a_dat, c_dat, e_dat, 'o', markersize=4, color='k') axs.plot_surface(A, C, E, linewidth=0, cmap='rainbow', alpha=0.5) plt.tight_layout()
plt.show()
Least squares fitting by Python
Linear/Nonlinear regression, multi‐variate (surface fitting)
scipy.optimize.curve_fit()
Machine learning (brief introduction) Machine learning:
Improve algorithm automatically through experience
Supervised learning
Input {x(0) }
Goal: to determine
Output {y(0) }
Experience: input‐output pairs
E.g. Least squares fitting
Function,
{x(1)} Operator, {y(1)}
Unsupervised learning
Recognize a pattern Only input
E.g. Clustering data
i
i
i Neural network i
Supervised learning
Experience: input‐output pairs E.g. Least squares fitting
Machine learning (brief introduction)
Usually, machine learning involves the problem of classification according to human’s experience, e.g. {Yes, No}, {True, False}, {1, 2, 3, ...}, {red, blue, green, ...}
E.g.
Input A lot of photos
Output If the photo shows a cat
A lot of photos of handwriting numbers
Recognize Arabic numbers
[Y. LeCun, L. Bottou, Y. Bengio, P. Haffner, Proceedings of the IEEE 86, 1998]
Coordinate of atoms of amorphous materials
If structural instability occurs locally
[E.D. Cubuk, S.S. Schoenholtz, J.M. Rieser, et. al. Phys. Rev. Lett. 114, 2015]
Input: {xi(0)}, {xi(1)}, ... ; Ouput: {yi}, yi = 0 or 1 (binary) Hypothesis (the function to be fitted):
g(z) 1 1ez
z
In neural network, ReLU function is used: g(z)ReLU(z)max(0,z)
Logistic regression
h(x;w)g w(0) w(1)x(1) w(N)x(N) g bw(1)x(1) w(N)x(N)
h(x;w,b) gwT x b g(x) can be logistic function:
g
1
O
g
O
z
h(x;w,b) gwT x b g(z) 1
g
parameters are w and b
h(x;w,b) Pr(y 1|x;w,b)
Logistic regression
1ez
Meaning of h(x; w, b): probability of y = 1 when the input is x and the
z
y1, h0.5
0, h 0.5
1
O
Logistic regression
In fitting problem, we need to construct cost function J = ? (e.g. least squares of residual) to minimize
Maximum Likelihood Estimation (MLE)
Probability that the outputs are y = (y1, ..., yM)T when the inputs are X = {xi(n)} and the parameters are w = (w(1), ..., w(N))T and b
Pr(y|X;w,b)
We hope that, for a certain (w, b), when the inputs are X, probably the outputs are y. MLE: determine the parameters (w, b) such that this probability is maximized
i i
i
i
J(w,b) 1 lnPr(y|X;w,b) M
J looks like entropy in form: S plnp
Logistic regression
The probability of y = 1 is h; the probability of y = 0 is 1 ‒ h
h(x ;w,b), y 1 Pr(y |x ;w,b) i i
y 1y
1h(x;w,b), y0 ii
h(x ;w,b) i 1h(x ;w,b) i
M
Pr(y|X;w,b) Pr(y |x ;w,b)
ii i1
M y ln h(x ;w,b) (1y )ln 1h(x ;w,b) 1M
iiii i1
ii i
Minimizing J is equivalent to maximizing cross entropy
Again, we can minimize J w.r.t. w and b to determine w and b. Compare least squares fitting and logistic regression: different J is used
wn1 wn α J
w wwn
w n
αM α inniin nn
Logistic regression
Similar to nonlinear regression, minimize J by gradient descent (iterative method)
bbn
h(x;w ,b)y x w XT h(X;w ,b)y
M i1 M b bαJ
n1 n
b wwn bbn
b n
h(x;w ,b)y b sumh(X;w ,b)y αM α
innin nn M i1 M
Receive signal
0: Not activated >0: Activated
x(1) x(2) x(3)
ReLU
Logisc regression → A neuron
z= g(z) wTx+b
h
Activated or not
Input
Hidden layers Layer [2]
Output Layer [3]
x(1)
a(1)[0]
w(1)[2] b(1)[2]
g a(1)[2]
x
a
x
a
(2)
(2)[0]
w(2)[2] b(2)[2]
g a(2)[2]
(3)
(3)[0]
w(3)[2] b(3)[2]
g a(3)[2]
w(1)[1] b(1)[1]
g a(1)[1]
w(2)[1] b(2)[1]
g a(2)[1]
w(1)[3] g a(1)[3] h(1) b(1)[3]
w(3)[1] b(3)[1]
g a(3)[1]
w(2)[3] g a(2)[3] h(2) b(2)[3]
w(4)[1] b(4)[1]
g a(4)[1]
Neural network (deep learning)
Layer [1]
Cost function for logistic regression: 1M
Cost function for neural network:
J()
y(k)lnh(k)(x ;) 1y(k) ln1h(k)(x ;)
1MK
b(i)[l] b(i)[l] α J
n1 n
b(i)[l]
Neural network
J(w,b)M y ln h(x ;w,b) (1y )ln 1h(x ;w,b)
iiii i1
Mi1k1
iiii w(ij)[l] ,b(i)[l]
l = 1, …, number of layers
For each l, i = 1, …, number of nodes of the l‐th layer
For each l, j = 1, …, number of nodes of the (l ‒ 1)‐th layer Gradient descent:
w(ij)[l] w(ij)[l] α J
Turing Award 2018: Yoshua Bengio, Geoffrey Hinton and Yann LeCun for making deep neural network a critical component of computing
n1 n
w(ij)[l]
n
n
Data types and methods
Recap
1. Basic programming I: Python
int, float, bool, str
Array of data: list, tuple, dictionary, set; related methods NumPy array
Variables and constants: evaluate and check types Looping: for, while; know how to write code block Condition: if, elif, else ; know how to write code block Branching: break, continue
Functions: define and call functions
File handling: open, read, write and close a file
NumPy
2. Basic programming II: NumPy, matplotlib, SymPy
matplotlib
Recap
Create different arrays
Check dimension, check/change type; index, search and sort arrays Element‐wise operations
Linear algebra
Vectors and matrices represented by NumPy arrays
inner product, cross product, matrix multiplication, etc.
Curve plot, polar plot, plot data from files, histogram, contour plot, 3D plot (curves and surfaces in 3D), subplots
Open, read and plot an image
SymPy
Symbols, expressions, substitution, convert to python functions, derivatives, integrals, limits, series expansion, summation of sequence, solve equations, solve differential equations, etc.
Gaussian elimination
Python code
Computational complexity
Recap
3. Systems of linear equations I: Direct method
Count the number of multiplications for input with size N LU decomposition
What is LU decomposition? Python code, SciPy function
QR decomposition
What is QR decomposition? Python code, SciPy function Compare the computation times for different algorithms How to create random matrices? (see Assignment 4 Q4)
Definitions and NumPy functions
Contraction
Priori and posteriori errors
Jacobi iteration: Python code Gauss‐Seidel iteration: Python code
Recap
4. Systems of linear equations II: Iterative method
When should we use iterative method? Vector norm and matrix norm
GMRES and conjugate gradient: SciPy functions
How to create sparse matrices? Convert NumPy array to compressed sparse row matrix (csr) or convert csr to NumPy array
Fixed‐point problem
Can fixed point be found by iteration? Schematic of procedure
SciPy function
Root‐finding problem
Netwon’s method: Schematic of procedure Damped Newton method: Schematic of procedure Chord method: Schematic of procedure
SciPy function
Recap
5. Systems of nonlinear equations
Linear regression
Line fitting
Nonlinear regression
Recap
6. Least squares fitting
Curve fitting (single‐variate), surface fitting (multi‐variate) Implementation by SciPy functions