|
利用海森矩阵判定多元函数的极值
海森矩阵（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))
|