141 凸最適化

Consider following quadratic term:

\begin{align} f(x) &= \frac{1}{2}x^TQx + bx, \\ \nabla f(x) &= Qx + b, \\ \nabla^2 f(x) &= Q. \end{align}

読み込み

In [1]:
import sympy as sy
sy.init_printing(use_latex=True)

from scipy.optimize import minimize
import numpy as np

import pandas as pd
import multiprocessing as mp

import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
% matplotlib inline

import warnings
warnings.filterwarnings("ignore")

変数定義

In [2]:
x_mat =  sy.Matrix(2, 1, sy.symbols('x:2:1'))
Q_mat = sy.Matrix(2, 2, sy.symbols('Q:2:2'))
b_mat = sy.Matrix(1, 2, sy.symbols('b:1:2'))
f_mat = sy.Function('f_mat')(x_mat)
f_mat = ((1/2) *x_mat.T*Q_mat*x_mat+b_mat*x_mat)
display(x_mat, Q_mat, b_mat)
display(Q_mat.det())
display(f_mat)
$$\left[\begin{matrix}x_{00}\\x_{10}\end{matrix}\right]$$
$$\left[\begin{matrix}Q_{00} & Q_{01}\\Q_{10} & Q_{11}\end{matrix}\right]$$
$$\left[\begin{matrix}b_{00} & b_{01}\end{matrix}\right]$$
$$Q_{00} Q_{11} - Q_{01} Q_{10}$$
$$\left[\begin{matrix}b_{00} x_{00} + b_{01} x_{10} + x_{00} \left(0.5 Q_{00} x_{00} + 0.5 Q_{10} x_{10}\right) + x_{10} \left(0.5 Q_{01} x_{00} + 0.5 Q_{11} x_{10}\right)\end{matrix}\right]$$

解析解へ具体的な値を代入

Optimum vector \(x^*\) is

\begin{align} \nabla f(x^*) &= Qx^* + b = 0, \\ x^* &= -Q^{-1}b. \end{align}
In [3]:
def isPSD(syA, tol=1e-8):
    eigenvals = syA.eigenvals()
    for key, value in eigenvals.items():
        if not key >= tol:
            yield False
    return True

Q_params = [(Q_mat[0,0], 1), (Q_mat[0,1], -1),
            (Q_mat[1,0], -1), (Q_mat[1,1], 2)
            ]
Q_mat_with_param = Q_mat.subs(Q_params)
if isPSD(Q_mat_with_param):
    Qinv = Q_mat_with_param.inv()
else:
    print('The matrix is not positive semidefinite symmetric matrix')

b_params = [(b_mat[0,0], 1), (b_mat[0,1], 1)]
b_mat_with_param = b_mat.subs(b_params)

x_optimal_analytical = -1*Qinv * b_mat_with_param.T
x_params = [(x_mat[0,0], x_optimal_analytical[0]),
            (x_mat[1,0], x_optimal_analytical[1])]
params = Q_params
params.extend(b_params)
params.extend(x_params)
f_optimal_analytical = f_mat.subs(params)[0].expand().simplify()

print('Q')
display(Q_mat_with_param)
print('Q inv')
display(Qinv)
print('b.T')
display(b_mat_with_param.T)
print('Optimal x')
display(x_optimal_analytical)
print('f(x*)')
display(f_optimal_analytical)
Q
$$\left[\begin{matrix}1 & -1\\-1 & 2\end{matrix}\right]$$
Q inv
$$\left[\begin{matrix}2 & 1\\1 & 1\end{matrix}\right]$$
b.T
$$\left[\begin{matrix}1\\1\end{matrix}\right]$$
Optimal x
$$\left[\begin{matrix}-3\\-2\end{matrix}\right]$$
f(x*)
$$-2.5$$

代数的にJacobian, Hessianを関数として作成

In [4]:
matrix_params = [(b_mat[0,0], 1), (b_mat[0,1], 1),
                 (Q_mat[0,0], 1), (Q_mat[0,1], -1),
                 (Q_mat[1,0], -1), (Q_mat[1,1], 2),
                ]
f_mat_with_param = f_mat.subs(matrix_params)[0].expand().simplify()

np_func = sy.lambdify(x_mat, f_mat_with_param, modules='numpy')

jac_f = [f_mat_with_param.diff(x) for x in x_mat]
jac_fn = [sy.lambdify(x_mat, jf, modules='numpy') for jf in jac_f]

hess_f = [[jac.diff(xx) for jac in jac_f] for xx in x_mat]
hess_fn = [sy.lambdify(x_mat, hf, modules='numpy') for hf in hess_f]

display(jac_f, jac_fn, hess_f, hess_fn)
$$\left [ 1.0 x_{00} - 1.0 x_{10} + 1, \quad - 1.0 x_{00} + 2.0 x_{10} + 1\right ]$$
[<function numpy.<lambda>>, <function numpy.<lambda>>]
$$\left [ \left [ 1.0, \quad -1.0\right ], \quad \left [ -1.0, \quad 2.0\right ]\right ]$$
[<function numpy.<lambda>>, <function numpy.<lambda>>]

Scipy.optimizeで数値的に最適化

jacobian,hessian関数をラップして変数を与える関数を作成

In [5]:
numerical_results = {}
init_guess = [0, 0]

def wrap_np_func(x):
    return np_func(*tuple(x))

def wrap_jacobian(x):
    return np.array([jfn(*tuple(x)) for jfn in jac_fn])

def wrap_hessian(x):
    return np.array([hfn(*tuple(x)) for hfn in hess_fn])

def cbf(x):
    global count
    count += 1
    f = wrap_np_func(x)
    plt.scatter(count, f)

Nelder-Mead

https://docs.scipy.org/doc/scipy/reference/optimize.minimize-neldermead.html#optimize-minimize-neldermead

In [6]:
method = 'Nelder-Mead'
count = 0
result = minimize(wrap_np_func, init_guess, callback=cbf,
                  method=method)

x_optimal_numerical = np.array(result.x)
x_params = [(x_mat[0,0], x_optimal_numerical[0]),
            (x_mat[1,0], x_optimal_numerical[1])]
numerical_results[method] = {'x':x_optimal_numerical,
                             'y': f_mat_with_param.subs(x_params)}

print('Optimal x')
display(numerical_results[method]['x'])
print('f(x*)')
display(numerical_results[method]['y'])
Optimal x
array([-2.99997629, -2.00000619])
f(x*)
$$-2.49999999953364$$
../_images/materials_141_Convex_optimization_13_4.png

Powell

https://docs.scipy.org/doc/scipy/reference/optimize.minimize-powell.html#optimize-minimize-powell

In [7]:
method = 'Powell'
count = 0
result = minimize(wrap_np_func, init_guess, callback=cbf,
                  method=method)

x_optimal_numerical = np.array(result.x)
x_params = [(x_mat[0,0], x_optimal_numerical[0]),
            (x_mat[1,0], x_optimal_numerical[1])]
numerical_results[method] = {'x':x_optimal_numerical,
                             'y': f_mat_with_param.subs(x_params)}

print('Optimal x')
display(numerical_results[method]['x'])
print('f(x*)')
display(numerical_results[method]['y'])
Optimal x
array([-3.        , -2.00000001])
f(x*)
$$-2.5$$
../_images/materials_141_Convex_optimization_15_4.png

Newton-CG

https://docs.scipy.org/doc/scipy/reference/optimize.minimize-newtoncg.html#optimize-minimize-newtoncg

In [8]:
method = 'Newton-CG'
count = 0
result = minimize(wrap_np_func, init_guess, callback=cbf,
                  method=method, jac=wrap_jacobian)

x_optimal_numerical = np.array(result.x)
x_params = [(x_mat[0,0], x_optimal_numerical[0]),
            (x_mat[1,0], x_optimal_numerical[1])]
numerical_results[method] = {'x':x_optimal_numerical,
                             'y': f_mat_with_param.subs(x_params)}

print('Optimal x')
display(numerical_results[method]['x'])
print('f(x*)')
display(numerical_results[method]['y'])
Optimal x
array([-3., -2.])
f(x*)
$$-2.5$$
../_images/materials_141_Convex_optimization_17_4.png

trust-krylov

https://docs.scipy.org/doc/scipy/reference/optimize.minimize-trustkrylov.html

In [9]:
method = 'trust-krylov'
count = 0
result = minimize(wrap_np_func, init_guess, callback=cbf,
                  method=method, jac=wrap_jacobian, hess=wrap_hessian)

x_optimal_numerical = np.array(result.x)
x_params = [(x_mat[0,0], x_optimal_numerical[0]),
            (x_mat[1,0], x_optimal_numerical[1])]
numerical_results[method] = {'x':x_optimal_numerical,
                             'y': f_mat_with_param.subs(x_params)}

print('Optimal x')
display(numerical_results[method]['x'])
print('f(x*)')
display(numerical_results[method]['y'])
Optimal x
array([-3., -2.])
f(x*)
$$-2.5$$
../_images/materials_141_Convex_optimization_19_4.png

解の可視化

関数へ数値の代入

In [10]:
def calc_func(arg):
    x1, x2 = arg
    return [x1, x2, np_func(x1, x2)]

def calc_parallel(func, x_lists, column_list):
    df = pd.DataFrame(columns=column_list)

    pool = mp.Pool(processes=mp.cpu_count() - 1)
    tmp_df = pd.DataFrame(pool.map(func, x_lists), columns=column_list)
    pool.close()
    pool.join()
    return df.append(tmp_df)

x_min = -10
x_max = 10
x_delta = 1
x_lists = [ (float(x1), float(x2))
                 for x1 in list(range(x_min, x_max, x_delta))
                 for x2 in list(range(x_min, x_max, x_delta))
             ]
column_list = ['x1', 'x2', 'y']

func_df = calc_parallel(calc_func, x_lists, column_list)

グラフ化

In [11]:
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

func_pivot_df = func_df.pivot(index='x1', columns='x2', values='y')
x1_pivot = func_pivot_df.index.values
x2_pivot = func_pivot_df.columns.values
y_pivot = func_pivot_df.values

ax.plot_wireframe(x1_pivot, x2_pivot, y_pivot)
ax.scatter(x_optimal_analytical[0], x_optimal_analytical[1], f_optimal_analytical,
           s=100, c='r', label='Analitic')
for key, value in numerical_results.items():
    ax.scatter(value['x'][0], value['x'][1], value['y'],
               s=100, label='{}'.format(key))

ax.set_xlabel('$x_1$')
ax.set_ylabel('$x_2$')
ax.set_zlabel('$f(x)$', rotation=90)
plt.legend(bbox_to_anchor=(1.05, 1), loc='upper left', borderaxespad=0)
plt.show()
../_images/materials_141_Convex_optimization_23_0.png