▼ Calculus (미적분)
SciPy
- 파이썬 기반의 오픈소스 계산 툴이다.
Matpolotlib.pyplot
- Equation
import matplotlib.pyplot as plt
xs = [x for x in range(-4, 10)]
ys = [0.1*x**3 - 0.8*x**2 - 1.5*x + 5.4 for x in xs]
plt.plot(xs, ys, 'r-')
plt.show()
- plot([x], y, [fmt], …)
: 이 함수는 x에 대한 y를 선과/또는 마커로 표시. 'fmt' 문자열을 통해 색상, 마커 스타일, 선 유형 등의 기본 형식을 정의 가능. - hist(x, bins=None, range=None, …)
: 이 함수는 히스토그램을 그려낸다.. 배열 x를 받아 데이터를 구간별로 나눕니다. 구간의 수, 범위 등 다양한 속성을 지정 가능. - contour([X, Y,] Z, [levels], …)
: 이 함수는 Z 데이터의 등고선을 그립니다. 3차원 데이터를 2차원 플롯에 그리는 데 유용하다.
( X와 Y는 격자를 정의 ! ) - imshow(X, cmap=None, …)
: 이 함수는 데이터를 2차원 정규 래스터(raster) 형태의 이미지로 표시하고, cmap을 통해 색상 맵을 설정 가능.
( raster: 많은 양의 픽셀(Pixel)이 모여 하나의 이미지를 구성하는 방식 ) - title(label, …)
: 이 함수는 축(axis)에 제목을 설정. - axis(…)
: 이 함수는 축 속성(예: 'on'/'off', 'equal', 축 범위 등)을 가져오거나 설정. - legend(…)
: 이 함수는 축에 legend(선들을 구분해주기 위한 이름)를 배치. - grid(…)
: 이 함수는 grid(격자선)을 설정. - xlabel(label, …), xlim(left, right), …
: 이 함수들은 x축의 label과 limit를 설정. - figure(num=None, …)
: 이 함수는 새로운 그림을 생성하거나 기존 그림을 활성화. - show(…)
: 이 함수는 모든 열린 그림을 표시. - savefig(filename, …)
: 이 함수는 현재 그림을 파일로 저장.
- Equation
import matplotlib.pyplot as plt
scale = 10
xs = [x/scale for x in range(-4*scale, 10*scale)]
ys = [0.1*x**3 - 0.8*x**2 - 1.5*x + 5.4 for x in xs]
plt.title('$y = 0.1x^3 - 0.8x^2 - 1.5x + 5.4$') # $ -> LaTeX style
plt.plot(xs, ys, 'r-')
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.axis('equal')
plt.show()
Differentiation (미분)
- 주어진 함수의 도함수(기울기) 혹은 선형 근사(linear approximation)를 찾는 과정
import matplotlib.pyplot as plt
scale = 10
xs = [x/scale for x in range(-4*scale, 10*scale)]
ys = [0.1*x**3 - 0.8*x**2 - 1.5*x + 5.4 for x in xs]
yt = [-3.5*x + 7 for x in xs]
plt.plot(xs, ys, 'r-', label='y')
plt.plot(xs, yt, 'b--', label='tangent line at x=2')
plt.xlabel('x')
plt.ylabel('y')
plt.grid()
plt.legend()
plt.show()
- SymPy - 도함수 구하기
import sympy as sp
x, y = sp.symbols('x y')
y = 0.1*x**3 - 0.8*x**2 - 1.5*x + 5.4
yd = sp.diff(y, x) // dy를 dx로 나누기 (미분)
print(yd) # 0.3*x**2 - 1.6*x - 1.5
print(float(yd.subs({x: 2}))) # -3.5 / x를 2로 치환하고 float타입으로 출력
- SymPy - y가 0일 때 x 찾기
roots = sp.solveset(y, x) # y: 함수, x: 구하고 싶은 미지수
print(roots) # FiniteSet(-3.0, 2.0, 9.0)
r0 = float(roots.args[0]) # -3.0 / Casting from Float object to float
- 인수 분해 (Factorize)
sp.pprint(sp.factor(y)) # 5.4·(0.11...·x - 1.0)·(0.33...·x + 1.0)·(0.5·x - 1.0)
y = (x**3 - 8*x**2 - 15*x + 54)/10 # 분수 형태로 변환
sp.pprint(sp.factor(y)) # (x - 9)·(x - 2)·(x + 3)
# ───────────────────────
# 10
print(sp.factor(y)) # (x - 9)*(x - 2)*(x + 3)/10
- 함수가 연속인 경우엔 아래처럼 기울기를 쉽게 구해볼 수 있다.
y = lambda x: 0.1*x**3 - 0.8*x**2 - 1.5*x + 5.4
x0 = 2
h = 0.001
d1 = (y(x0+h) - y(x0)) / h # -3.5002
d2 = (y(x0+h/2) - y(x0-h/2)) / h # -3.5000
Integration (적분)
- 잘 사용 안한다 !
▼ Linear Algebra
Vector
- (Physical quantities) 속도와 방향을 모두 표현하기 위해 사용 !
- A point [1, 4] , a line 2x + 3y − 14 = 0 → [2, 3, −14] , and its normal vector [2, 3] or [0.554 … , 8.832 …] 로 표현하기 위해 !
Matrix
- 관성 모멘트(Moment of inertia) 같은 physical quantities를 표현하기 위해 사용 !
- 데이터와 모델(model)을 표현하기 위해 !
ex) Images, multiple vectors - geometric tranformation을 표현하기 위해 !
- System of linear equations(선형연립방정식)을 풀기 위해서 !
- Matrix multiplication
➜ Not commutative (교환법칙 성립 X)
import sympy as sp
a11, a12, a21, a22 = sp.symbols('a11 a12 a21 a22')
b11, b12, b21, b22 = sp.symbols('b11 b12 b21 b22')
A = sp.Matrix([[a11, a12], [a21, a22]])
B = sp.Matrix([[b11, b12], [b21, b22]])
print(A+B == B+A) # True
print(A*B == B*A) # False
print(A*B - B*A) # Check their difference
- Matrix inverse (AA^-1 = A-1A = I; 역수)
- Pseudo-inverse (기존 Matrix inverse를 확장한 개념)
➜ 반드시 square matrix일 필요 X
NumPy
- 다차원 배열과 행렬을 계산할 수 있는 함수를 제공하는 중요한 라이브러리이다 !
- numpy.array
➜ homogenous(~same) data type을 가진다. (vs. list and tuple)
import numpy as np
# 1. Create an array from a composite data (list or tuple, not set and dictionary)
A = np.array([3, 29, 82])
B = np.array(((3., 29, 82), (10, 18, 84)))
C = np.array([[3, 29, 82], [10, 18, 84]], dtype=float)
D = np.array([3, 29, 'Shin'])
E = np.array([[3], [29], [82]])
# nidm: 몇차원, size: 크기, shape: n x m, dtype: 데이터 타입
print(A.ndim, A.size, A.shape, A.dtype) # 1 3 (3,) int32
# Note) np.array_equal(A, A.T) == True
# 1차원 배열은 전치시켜도 같다.
print(B.ndim, B.size, B.shape, B.dtype) # 2 6 (2, 3) float64
print(C.ndim, C.size, C.shape, C.dtype) # 2 6 (2, 3) float64
print(D.ndim, D.size, D.shape, D.dtype) # 1 3 (3,) <U11 (Unicode)
# Note) array(['3', '29', 'Shin'])
print(E.ndim, E.size, E.shape, E.dtype) # 2 3 (3,1) int32
# Note) np.array_equal(E, E.T) == False
# because E.T == array([[3, 29, 82]])
# 2. Create an array using initializers
F = np.zeros((3, 2)) # Create a 3x2 array filled with 0 (default: float64)
G = np.ones((3, 2)) # Create a 2x3 array filled with 1
H = np.eye(3, dtype=np.float32) # Create a 3x3 identity matrix (single-precision)
I = np.empty((3, 2)) # == np.zeros((3, 2))
J = np.empty((0, 9)) # [ ] with size of (0, 9)
K = np.arange(0, 1, 0.2) # 간격으로 쪼갠다 - Step 0.2: array([0., 0.2, 0.4, 0.6, 0.8])
L = np.linspace(0, 1, 5) # 개수로 쪼갠다 - Number 5: array([0., 0.25, 0.50, 0.75, 1.])
M = np.random.random((3, 2)) # == np.random.uniform(size=(3, 2))
# Note) np.random.normal()
- Indexing & Slicing
import numpy as np
A = np.array(((3., 29, 82), (10, 18, 84)))
# 1. Indexing and slicing
A[1][1] # 18.0
A[1, 1] # 18.0 Note) list[1, 1] does not work!
A[1,1:2] # array([18.])
A[1,:] # Get a row: array([10., 18., 84.])
A[:,2] # Get a column: array([82., 84.])
A[0:2,0:2] # Get a submatrix: array([[3., 29.], [10., 18.]])
# 2. Logical indexing
A > 80 # array([[False, False, True], [False, False, True]])
A[A > 80] # array([82., 84.])
A[A > 80] = 80 # Masked operations are possible!
# 3. Fancy indexing
A[(1, 0, 0), (1, 0, 2)] # array([18., 3., 82.])
# Get items at (1, 1), (0, 0), (0, 2)
- Arithmetic operations (cf. different form list)
import numpy as np
A = np.array([[1, 2], [3, 4]])
B = np.array([[5, 6], [7, 8]])
# 1. Element-wise arithmetic operations
print(A + B) # np.add(A, B)
print(A - B) # np.subtract(A, B)
print(A * B) # np.multiply(A, B)
print(A / B) # np.divide(A, B)
# 2. Matrix operations
print(A.T) # A.transpose()
print(A @ B) # np.matmul(A, B) : Matrix multiplication
print(np.linalg.norm(A)) # 5.48 (default: L2-norm)
print(np.linalg.matrix_rank(A)) # 2, full rank
print(np.linalg.det(A)) # -2, non-zero determinant
print(np.linalg.inv(A)) # Matrix inverse
print(np.linalg.pinv(A)) # Matrix pseudo-inverse Note) inv()의 기능을 포함하고 있다.
# 3. Broadcasting
print(A + 1) # [[2, 3], [4, 5]]
print(A + [0, -1]) # [[1, 1], [3, 3]]
print(A + [[1], [-1]]) # [[2, 3], [2, 3]]
- Broatcasting 개념
- Example - Line fitting from two points, (1, 4) and (4, 2)
import numpy as np
A = np.array([[1., 1.], [4., 1.]])
b = np.array([[4.], [2.]])
A_inv = np.linalg.inv(A)
print(A_inv * b) # [[-1.33333333 1.33333333]
# [ 2.66666667 -0.66666667]] Note) broadcast
print(A_inv @ b) # [[-0.66666667]
# [ 4.66666667]]
- Example - Line fitting from more than two points such as (1, 4), (4, 2) and (7,1)
import numpy as np
A = np.array([[1., 1.], [4., 1.], [7., 1.]])
b = np.array([[4.], [2.], [1.]])
A_inv = np.linalg.pinv(A) # Left inverse
print(A_inv @ b) # [[-0.5]
# [ 4.33333333]]
- Example - Line fitting with noisy but more data
➜ 데이터가 많을수록 'noise' 에 강해진다.(robust)
➜ 'Line fitting' 이 잘되는 것을 확인 가능 !
- Curve fitting (≒ Line fitting)
import numpy as np
import matplotlib.pyplot as plt
true_curve = lambda x: 0.1*x**3 - 0.8*x**2 - 1.5*x + 5.4
data_range = (-6, 12)
data_num = 100
noise_std = 0.5
# Generate the true data
x = np.random.uniform(data_range[0], data_range[1], size=data_num)
y = true_curve(x)
# Add Gaussian noise
xn = x + np.random.normal(scale=noise_std, size=x.shape)
yn = y + np.random.normal(scale=noise_std, size=y.shape)
# Solve the system of equations
A = np.vstack((xn**3, xn**2, xn, np.ones(xn.shape))).T
b = yn
curve = np.linalg.pinv(A) @ b
# Plot the data and result
plt.title(f'Curve: y={curve[0]:.3f}*$x^3$ + {curve[1]:.3f}*$x^2$ + {curve[2]:.3f}*$x$ + {curve[3]:.3f}')
xc = np.linspace(*data_range, 100)
plt.plot(xc, true_curve(xc), 'r-', label='The true curve')
plt.plot(xn, yn, 'b.', label='Noisy data')
plt.plot(xc, curve[0]*xc**3 + curve[1]*xc**2 + curve[2]*xc + curve[3], 'g-', label='Estimate')
plt.xlim(data_range)
plt.legend()
plt.show()
- Curve fitting + Model selection
import numpy as np
import matplotlib.pyplot as plt
# 모델 A를 만들어주는 기능을 추가해준 것 뿐이다.
def buildA(order, xs):
A = np.empty((0, len(xs)))
for i in range(order + 1):
A = np.vstack((xs**i, A))
return A.T
true_coeff = [0.1, -0.8, -1.5, 5.4]
poly_order = 3 # Try other integer (>= 0)
data_range = (-6, 12)
data_num = 100
noise_std = 1
# Generate the true data
x = np.random.uniform(data_range[0], data_range[1], size=data_num)
y = buildA(len(true_coeff) - 1, x) @ true_coeff
# Add Gaussian noise
xn = x + np.random.normal(scale=noise_std, size=x.shape)
yn = y + np.random.normal(scale=noise_std, size=y.shape)
# Solve the system of equations
A = buildA(poly_order, xn)
b = yn
coeff = np.linalg.pinv(A) @ b
# Plot the data and result
plt.title(f'Order: {poly_order}, Coeff: ' + np.array2string(coeff, precision=2, suppress_small=True))
xc = np.linspace(*data_range, 100)
plt.plot(xc, np.matmul(buildA(len(true_coeff) - 1, xc), true_coeff), 'k-', label='The true curve', alpha=0.2)
plt.plot(xn, yn, 'b.', label='Noisy data')
plt.plot(xc, np.matmul(buildA(poly_order, xc), coeff), 'g-', label='Estimate')
plt.xlim(data_range)
plt.legend()
plt.show()
- inverse가 안되는 경우
➜ 'inv', 'pinv' 모두 오류 발생 !
(해결 방법)
➜ 아래에서 v로 구성된 것을 "Null space" 라 하는데, 이 'Null' 벡터를 찾으면 된다.
import numpy as np
A = np.array([[1., 1.], [1., 1.]])
b = np.array([[4.], [2.]])
A_inv = np.linalg.inv(A) # Error! (singular matrix)
print(A_inv @ b)
# Null space를 찾는 방법
import numpy as np
from scipy import linalg
A = np.array([[1., 4., 1.], [1., 2., 1.]])
x = linalg.null_space(A)
print(x / x[0]) # [[1.], [0.], [-1.]] Note) Line: x – 1 = 0
▼ Optimization
Geometric distance
Optimization
- 'Objective functioin' 을 최대/최소화(minimization)시켜서 지표(criterion)상 가장 적합한(best) 것을 선택하는 것.
➜ 예를 들어 loss function을 어떤 것을 선택하는지에 따라 성능이 매우 달라질 수 있다 !
Nonlinear Optimization
- Gradient descent
➜ 기울기가 감소하는 방향으로 !
➜ 미분 가능한 함수의 local minimum(극소값)을 찾기 위해 1차 편미분을 반복적으로 하며 update하는 방법.
import numpy as np
import matplotlib.pyplot as plt
f = lambda x: 0.1*x**3 - 0.8*x**2 - 1.5*x + 5.4
fd = lambda x: 0.3*x**2 - 1.6*x - 1.5
viz_range = np.array([-6, 12])
learn_rate = 0.1 # step size, Try 0.001, 0.01, 0.5, and 0.6 :
max_iter = 100
min_tol = 1e-6
x_init = 12 # 시작점, Try -1 -> 음의 무한대로 간다.
# Prepare visualization
xs = np.linspace(*viz_range, 100)
plt.plot(xs, f(xs), 'r-', label='f(x)', linewidth=2)
plt.plot(x_init, f(x_init), 'b.', label='Each step', markersize=12)
plt.axis((*viz_range, *f(viz_range))) // axis : xlim과 ylim을 한 번에 설정
plt.legend()
x = x_init
for i in range(max_iter):
# Run the gradient descent
xp = x
x = x - learn_rate*fd(x) # x 값을 조금씩 감소시키면서 극소(local minimum)을 찾는다
# Update visualization for each iteration
print(f'Iter: {i}, x = {xp:.3f} to {x:.3f}, f(x) = {f(xp):.3f} to {f(x):.3f} (f\'(x) = {fd(xp):.3f})')
lcolor = np.random.rand(3)
approx = fd(xp)*(xs-xp) + f(xp)
plt.plot(xs, approx, '-', linewidth=1, color=lcolor, alpha=0.5)
plt.plot(x, f(x), '.', color=lcolor, markersize=12)
# Check the terminal condition
if abs(x - xp) < min_tol:
break
plt.show()
- Newton's method
➜ 미분 가능한 함수의 극값(최대값 또는 최소값) 또는 해를 찾는데 사용되는 기법.
➜ 'step size' 필요 X, 'step size' 가 이미 정해져있다.
➜ 1/fdd(x)
import numpy as np
import matplotlib.pyplot as plt
f = lambda x: 0.1*x**3 - 0.8*x**2 - 1.5*x + 5.4
fd = lambda x: 0.3*x**2 - 1.6*x - 1.5
fdd = lambda x: 0.6*x - 1.6
viz_range = np.array([-6, 12])
max_iter = 100
min_tol = 1e-6
x_init = 12 # Try -2, 0, and 16/6 (a saddle point)
# Prepare visualization
xs = np.linspace(*viz_range, 100)
plt.plot(xs, f(xs), 'r-', label='f(x)', linewidth=2)
plt.plot(x_init, f(x_init), 'b.', label='Each step', markersize=12)
plt.axis((*viz_range, *f(viz_range)))
plt.legend()
x = x_init
for i in range(max_iter):
# Run the Newton method
xp = x
x = x - fd(x) / fdd(x)
# the maxima and saddle point problems를 해결하기 위한 방법
# ➜ Replace the denominator as abs(fdd(x)) and (abs(fdd(x)) + 1)
# Update visualization for each iteration
print(f'Iter: {i}, x = {xp:.3f} to {x:.3f}, f(x) = {f(xp):.3f} to {f(x):.3f} (f\'(x) = {fd(xp):.3f}, f\'\'(x) = {fdd(xp):.3f})')
lcolor = np.random.rand(3)
approx = 0.5*fdd(xp)*(xs-xp)**2 + fd(xp)*(xs-xp) + f(xp)
plt.plot(xs, approx, '-', linewidth=1, color=lcolor, alpha=0.8)
plt.plot(x, f(x), '.', color=lcolor, markersize=12)
# Check the terminal condition
if abs(x - xp) < min_tol:
break
plt.show()
scipy.optimize
- Optimization 문제나 Root finding에 대한 함수를 제공.
Minimization(최소화) ➜ 도함수를 구해줄 필요가 없다! ( 함수가 알아서 해주기 때문 )
➜ 파라미터 값을 조정해줄 수 있다. (tol, options)
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
f = lambda x: 0.1*x**3 - 0.8*x**2 - 1.5*x + 5.4
viz_range = np.array([-6, 12])
max_iter = 100
min_tol = 1e-6
x_init = 12 # Try -2, 0, and 16/6
# Find the minimum by SciPy
result = minimize(f, x_init, tol=min_tol, options={'maxiter': max_iter, 'return_all': True})
print(result)
# Visualize all iterations
xs = np.linspace(*viz_range, 100)
plt.plot(xs, f(xs), 'r-', label='f(x)', linewidth=2)
xr = np.vstack(result.allvecs)
plt.plot(xr, f(xr), 'b.', label='Each step', markersize=12)
plt.legend()
plt.axis((*viz_range, *f(viz_range)))
plt.show()
Line fitting with minimizing Algebraic and Geometric distance
- Algebraic distance
- Geometric distance는 선이 좀 더 가파른 기울기를 가질 때일수록 더 유용하게 쓰인다.
import numpy as np
import matplotlib.pyplot as plt
from scipy.optimize import minimize
true_line = lambda x: -14/3*x + 14/3
data_range = np.array([-4, 12])
data_num = 100
noise_std = 1
# Generate the true data
x = np.random.uniform(data_range[0], data_range[1], size=data_num)
y = true_line(x)
# Add Gaussian noise
xn = x + np.random.normal(scale=noise_std, size=x.shape)
yn = y + np.random.normal(scale=noise_std, size=y.shape)
# Find a line minimizing algebraic distance
A = np.vstack((xn, np.ones(xn.shape))).T
b = yn
l_alg = np.linalg.pinv(A) @ b
e_alg = np.mean(np.abs(l_alg[0]*xn - yn + l_alg[1]) / np.sqrt(l_alg[0]**2 + 1))
# Find a line minimizing geometric distance
geo_dist2 = lambda x: np.sum((x[0]*xn - yn + x[1])**2) / (x[0]**2 + 1)
result = minimize(geo_dist2, [-1, 0]) # The initial value: y = -x
l_geo = result.x
e_geo = np.mean(np.abs(l_geo[0]*xn - yn + l_geo[1]) / np.sqrt(l_geo[0]**2 + 1))
# Plot the data and result
plt.plot(data_range, true_line(data_range), 'r-', label='The true line')
plt.plot(xn, yn, 'b.', label='Noisy data')
plt.plot(data_range, l_alg[0]*data_range + l_alg[1], 'g-', label=f'Solve Ax=b (GeoError: {e_alg:.3f})')
plt.plot(data_range, l_geo[0]*data_range + l_geo[1], 'm-', label=f'Optimization (GeoError: {e_geo:.3f})')
plt.legend()
plt.xlim(data_range)
plt.show()
▼ Probablity
- 강의 자료를 참고하는 것이 좋을 것 같다.
