第十一篇： 利用海森矩陣判定多元函數的極值與最小二乘法

利用海森矩陣判定多元函數的極值

海森矩陣（Hessian Matrix），又譯作黑塞矩陣、海瑟矩陣、 海塞矩陣等，是一個多元函數的二階偏導數構成的方陣，描述 了函數的局部曲率。黑塞矩陣最早於19世紀由德國數學家 Ludwig Otto Hesse提出，並以其名字命名。海森矩陣常用於 解決優化問題，利用黑塞矩陣可判定多元函數的極值問題。

from sympy import pprint,Symbol,diff,solve,symarray,Eq,Expr,roots,simplify,lambdify,hessian
from sympy.matrices import Matrix,zeros,diag,eye,GramSchmidt
from sympy.abc import lamda,x,y,z,t
import numpy as np

def solve_extremum(expr,symList,stagnation_point:tuple=None):
    from sympy import hessian
    import numpy as np
    hsm = np.array(hessian(expr,symList).evalf()).astype(np.float32)
    sym_tup = tuple(symList)
    # 用特征值判斷
    eig_val = np.linalg.eigvals(hsm)
    print('\n海森矩陣的特征值：')
    print(eig_val)
    greater = eig_val > 0.0
    less = eig_val < 0.0
    bool_list = [True for i in range(len(eig_val))]
    if np.allclose(bool_list,greater):
        print('\nf{} = {}的海森矩陣是正定矩陣，存在極小值點。\n'.format(sym_tup, expr))
        if stagnation_point:
            print('故{}是極小值點，且極小值f{} = {}\n'.format(stagnation_point, stagnation_point,
                                                   expr.subs({k: v for k, v in zip(symList, stagnation_point)})))
    elif np.allclose(bool_list,less):
        print('\n' + 'f{} = {}的海森矩陣是負定矩陣，存在極大值點。\n'.format(sym_tup, expr))
        if stagnation_point:
            print('故{}是極大值點，且極大值f{} = {}\n'.format(stagnation_point, stagnation_point,
                                                   expr.subs({k: v for k, v in zip(symList, stagnation_point)})))
    else:
        print('\n無法判斷極值點\n')


'''三元函數'''
f = x**2+y**2+z**2+2*x+4*y-6*z

#求駐點
dx = diff(f,x)
dy = diff(f,y)
dz = diff(f,z)
stagnation_point = solve((dx,dy,dz),x,y,z)
if isinstance(stagnation_point,dict):
    point = []
    point.append((stagnation_point[x],stagnation_point[y],stagnation_point[z]))
    stagnation_point = point
for i,point in enumerate(stagnation_point):
    print('駐點{}：{}'.format(i+1,point))

f_hsm = hessian(f,[x,y,z]) #海森矩陣
print('\n'+'f(x,y,z)的海森矩陣：')
pprint(f_hsm)
solve_extremum(f,[x,y,z],stagnation_point[0])


A = Matrix([[-5,-2,4],[-2,-1,2],[4,2,-5]])
# A = Matrix([[-10,-4,-12],[-4,-2,14],[-12,14,-1]])
X = Matrix([x,y,z])
f2 = (X.T*A*X)[0,0].simplify()

#求駐點
dx = diff(f2,x)
dy = diff(f2,y)
dz = diff(f2,z)

stagnation_point = solve((dx,dy,dz),x,y,z)
if isinstance(stagnation_point,dict):
    point = []
    point.append((stagnation_point[x],stagnation_point[y],stagnation_point[z]))
    stagnation_point = point
for i,point in enumerate(stagnation_point):
    print('駐點{}：{}'.format(i+1,point))

f2_hsm = hessian(f2,[x,y,z,t])
print('\n'+'f(x,y,z)的海森矩陣：')
pprint(f2_hsm)
solve_extremum(f2,[x,y,z],stagnation_point[0])

'''隨機生成整數測試'''
# np.random.seed(28)
#生成在駐點附近的點
arr_x = np.random.uniform(-2,0,(5,))
arr_y = np.random.uniform(-3,-1,(5,))
arr_z = np.random.uniform(2,4,(5,))

#sympy中使用numpy高效計算
f = lambdify((x,y,z),f,'numpy')
print(f(arr_x,arr_y,arr_z))

f2 = lambdify((x,y,z),f2,'numpy')
print(f2(arr_x,arr_y,arr_z))

最小二乘法

最小二乘法

數據擬合

from sympy import pprint,Symbol,diff,solve,symarray,Eq,Expr,roots,simplify,lambdify
from sympy.matrices import Matrix,zeros,diag,eye,GramSchmidt
from sympy.abc import lamda,x,y,z,t
import numpy as np
A = Matrix([[2,1],[4,1],[6,1],[8,1],[10,1],[12,1]])
B = Matrix([8.9,10.1,11.2,12.0,13.1,13.9])
pprint(A.solve_least_squares(B))