# Optimization: principles and algorithms¶

Bierlaire, M. (2015). Optimization: Principles and Algorithms. EPFL Press.
The numbering of the algorithms, tables and page refer to the book.

## Chapter 14: Least squares problem¶

### Algorithm 14.1: Gauss-Newton¶

Solve $\min_x f(x)=\frac{1}{2} g(x)^T g(x)$ using $g(x)$ and $\nabla g(x)$.

In [1]:
import numpy as np
from scipy import linalg
# Define a function simple to call for the inner product
def scalarInner(a,b):
return(np.asscalar(a.T.dot(b)))

In [2]:
def gaussNewton(obj,x0,eps=1.0e-4,maxiter=100):
n = x0.shape[0]
xk = x0
iters = list()
iters.append([xk,0.5*scalarInner(g,g),linalg.norm(deriv)])
k = 0
while linalg.norm(deriv) > eps and k < maxiter:
xk = xk + dk ;
iters.append([xk,0.5*scalarInner(g,g),linalg.norm(deriv)])
k = k + 1
return xk,iters



Example 14.1, page 329

In [3]:
men = [ 77005.0,
76315.0,
70891.0,
67667.0,
64643.0,
61770.0,
61593.0,
63227.0,
63684.0,
66914.0,
72407.0,
82413.0,
86515.0,
84896.0,
79660.0,
75827.0,
72606.0,
69423.0,
69690.0,
69744.0,
70418.0,
71998.0,
77268.0,
87299.0]
women = [ 57312.0 ,
56839.0 ,
55501.0 ,
55491.0 ,
54217.0 ,
53098.0 ,
54701.0 ,
56596.0 ,
56663.0 ,
58622.0 ,
59660.0 ,
59896.0 ,
61643.0  ,
61105.0  ,
59333.0  ,
60024.0  ,
58684.0  ,
57075.0  ,
58826.0  ,
60212.0  ,
60654.0  ,
61445.0  ,
61805.0  ,
62138.0  ]

In [4]:
A = np.empty((len(men),2))
A[:,0] = np.array(men)
A[:,1] = 1.0
b = np.array(women).reshape(len(women),1)

def ex1401(x):
g = A.dot(x) - b
return g, A.T

In [5]:
x0 = np.array([[0],[0]])
(sol,iters) = gaussNewton(ex1401,x0)
sol

Out[5]:
array([[2.53435732e-01],
[3.99825112e+04]])
In [6]:
print("Number of iteration(s): {}".format(len(iters)-1))

Number of iteration(s): 1


Figure 14.1, page 331

In [7]:
%matplotlib inline
import matplotlib.pyplot as plt
t = np.arange(min(men),max(men), 100)
plt.plot(t,sol.item(0)*t + sol.item(1),color='red')
plt.scatter(men,women)
plt.show()


### Algorithm 14.2: Kalman filter¶

In [8]:
def kalmanFilter(x,H,A,b,lbd=1.0):
H = lbd * H + A.T.dot(A)
d = linalg.solve(H,A.T.dot(b-A.dot(x)))
x = x + d
return x,H


Example 14.5, page 340. Table 14.3.

In [9]:
A1 = A[0:12,:]
b1 = b[0:12]
x = np.zeros((2,1))
H = np.zeros((2,2))
x,H = kalmanFilter(x,H,A1,b1)
print("{}\t{}".format(x.item(0),x.item(1)))
for k in range(12,24):
x,H = kalmanFilter(x,H,A[k,:].reshape(1,2),b[k])
print("{}\t{}".format(x.item(0),x.item(1)))

0.21075800939512573	41998.073102822156
0.23963015306480753	40074.374324535864
0.24903841146458075	39451.46648857908
0.24935508197492404	39431.5502333248
0.25509174102240634	39122.56968410651
0.25545326724564976	39157.64458026606
0.25499801776510184	39200.62732569789
0.25108663539694964	39579.67951387378
0.24500964592951335	40172.30449301707
0.2409503118980672	40617.41272731181
0.2416309305444442	40726.71916388981
0.25301161807435313	40011.71877391723
0.2534357324827362	39982.51120367221