SciPy (phát âm là sai pie) là gói ứng dụng toán học dựa trên phần mở rộng Numpy Python. Với SciPy, phiên Python tương tác của bạn sẽ trở thành môi trường tạo nguyên mẫu hệ thống phức tạp và khoa học dữ liệu hoàn chỉnh giống như MATLAB, IDL, Octave, R-Lab và SciLab. Hôm nay tôi muốn nói ngắn gọn về cách sử dụng một số thuật toán tối ưu hóa nổi tiếng trong gói scipy.optimize. Bạn luôn có thể nhận được trợ giúp chi tiết và cập nhật hơn về cách sử dụng các hàm bằng cách sử dụng lệnh help() hoặc sử dụng Shift+Tab.
Giới thiệu
Để tránh cho bản thân và người đọc khỏi việc tìm kiếm và đọc các nguồn chính, các liên kết đến mô tả phương pháp sẽ chủ yếu có trên Wikipedia. Theo quy định, thông tin này đủ để hiểu các phương pháp nói chung và các điều kiện áp dụng chúng. Để hiểu bản chất của các phương pháp toán học, hãy theo các liên kết đến các ấn phẩm có thẩm quyền hơn, có thể tìm thấy ở cuối mỗi bài viết hoặc trong công cụ tìm kiếm yêu thích của bạn.
Vì vậy, mô-đun scipy.optizes bao gồm việc thực hiện các quy trình sau:
- Giảm thiểu có điều kiện và vô điều kiện các hàm vô hướng của một số biến (tối thiểu) bằng các thuật toán khác nhau (Nelder-Mead simplex, BFGS, gradient liên hợp Newton, и )
- Tối ưu hóa toàn cầu (ví dụ: , )
- Giảm thiểu dư lượng (least_squares) và thuật toán điều chỉnh đường cong sử dụng bình phương tối thiểu phi tuyến (curve_fit)
- Giảm thiểu hàm vô hướng của một biến (minim_scalar) và tìm kiếm gốc (root_scalar)
- Giải đa chiều hệ phương trình (gốc) bằng nhiều thuật toán khác nhau (hybrid Powell, hoặc các phương pháp quy mô lớn như ).
Trong bài viết này, chúng tôi sẽ chỉ xem xét mục đầu tiên trong toàn bộ danh sách này.
Giảm thiểu vô điều kiện hàm vô hướng của một số biến
Hàm tối thiểu từ gói scipy.optizes cung cấp giao diện chung để giải quyết các vấn đề giảm thiểu có điều kiện và vô điều kiện của hàm vô hướng của một số biến. Để chứng minh cách thức hoạt động của nó, chúng ta sẽ cần một hàm số phù hợp với nhiều biến mà chúng ta sẽ giảm thiểu theo nhiều cách khác nhau.
Với mục đích này, hàm Rosenbrock của N biến là hoàn hảo, có dạng:
Mặc dù thực tế là hàm Rosenbrock và các ma trận Jacobi và Hessian của nó (tương ứng là đạo hàm thứ nhất và thứ hai) đã được xác định trong gói scipy.optizes, chúng tôi sẽ tự xác định nó.
import numpy as np
def rosen(x):
"""The Rosenbrock function"""
return np.sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0, axis=0)Để rõ ràng, hãy vẽ 3D các giá trị của hàm Rosenbrock của hai biến.
Vẽ mã
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
# Настраиваем 3D график
fig = plt.figure(figsize=[15, 10])
ax = fig.gca(projection='3d')
# Задаем угол обзора
ax.view_init(45, 30)
# Создаем данные для графика
X = np.arange(-2, 2, 0.1)
Y = np.arange(-1, 3, 0.1)
X, Y = np.meshgrid(X, Y)
Z = rosen(np.array([X,Y]))
# Рисуем поверхность
surf = ax.plot_surface(X, Y, Z, cmap=cm.coolwarm)
plt.show()
Biết trước cực tiểu bằng 0 tại , hãy xem các ví dụ về cách xác định giá trị tối thiểu của hàm Rosenbrock bằng cách sử dụng các quy trình scipy.optizes khác nhau.
Phương pháp đơn giản Nelder-Mead
Giả sử có một điểm ban đầu x0 trong không gian 5 chiều. Hãy tìm điểm tối thiểu của hàm Rosenbrock gần nó nhất bằng thuật toán (thuật toán được chỉ định là giá trị của tham số phương thức):
from scipy.optimize import minimize
x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
res = minimize(rosen, x0, method='nelder-mead',
options={'xtol': 1e-8, 'disp': True})
print(res.x)Optimization terminated successfully.
Current function value: 0.000000
Iterations: 339
Function evaluations: 571
[1. 1. 1. 1. 1.]Phương pháp đơn hình là cách đơn giản nhất để giảm thiểu một hàm được xác định rõ ràng và khá trơn tru. Nó không yêu cầu tính đạo hàm của một hàm; chỉ cần xác định các giá trị của nó là đủ. Phương pháp Nelder-Mead là một lựa chọn tốt cho các bài toán tối thiểu hóa đơn giản. Tuy nhiên, vì nó không sử dụng ước tính độ dốc nên có thể mất nhiều thời gian hơn để tìm mức tối thiểu.
phương pháp Powell
Một thuật toán tối ưu hóa khác trong đó chỉ tính các giá trị hàm là . Để sử dụng nó, bạn cần đặt phương thức = 'powell' trong hàm minim.
x0 = np.array([1.3, 0.7, 0.8, 1.9, 1.2])
res = minimize(rosen, x0, method='powell',
options={'xtol': 1e-8, 'disp': True})
print(res.x)Optimization terminated successfully.
Current function value: 0.000000
Iterations: 19
Function evaluations: 1622
[1. 1. 1. 1. 1.]Thuật toán Broyden-Fletcher-Goldfarb-Shanno (BFGS)
Để đạt được sự hội tụ nhanh hơn về lời giải, thủ tục sử dụng độ dốc của hàm mục tiêu. Độ dốc có thể được chỉ định dưới dạng hàm hoặc được tính toán bằng cách sử dụng sai phân bậc nhất. Trong mọi trường hợp, phương thức BFGS thường yêu cầu ít lệnh gọi hàm hơn phương thức đơn giản.
Chúng ta hãy tìm đạo hàm của hàm Rosenbrock ở dạng phân tích:
Biểu thức này hợp lệ đối với đạo hàm của tất cả các biến ngoại trừ biến đầu tiên và biến cuối cùng, được định nghĩa là:
Hãy xem hàm Python tính toán độ dốc này:
def rosen_der (x):
xm = x [1: -1]
xm_m1 = x [: - 2]
xm_p1 = x [2:]
der = np.zeros_like (x)
der [1: -1] = 200 * (xm-xm_m1 ** 2) - 400 * (xm_p1 - xm ** 2) * xm - 2 * (1-xm)
der [0] = -400 * x [0] * (x [1] -x [0] ** 2) - 2 * (1-x [0])
der [-1] = 200 * (x [-1] -x [-2] ** 2)
return derHàm tính toán gradient được chỉ định là giá trị của tham số jac của hàm minim, như minh họa bên dưới.
res = minimize(rosen, x0, method='BFGS', jac=rosen_der, options={'disp': True})
print(res.x)Optimization terminated successfully.
Current function value: 0.000000
Iterations: 25
Function evaluations: 30
Gradient evaluations: 30
[1.00000004 1.0000001 1.00000021 1.00000044 1.00000092]Thuật toán gradient liên hợp (Newton)
Thuật toán là một phương pháp Newton đã được sửa đổi.
Phương pháp của Newton dựa trên việc xấp xỉ một hàm trong một khu vực cục bộ bằng đa thức bậc hai:
đâu là ma trận đạo hàm bậc hai (ma trận Hessian, Hessian).
Nếu Hessian là xác định dương thì có thể tìm thấy cực tiểu cục bộ của hàm này bằng cách đánh đồng gradient bằng XNUMX của dạng bậc hai bằng XNUMX. Kết quả sẽ là biểu thức:
Hessian nghịch đảo được tính bằng phương pháp gradient liên hợp. Một ví dụ về việc sử dụng phương pháp này để cực tiểu hóa hàm Rosenbrock được đưa ra dưới đây. Để sử dụng phương pháp Newton-CG, bạn phải chỉ định một hàm tính toán Hessian.
Hessian của hàm Rosenbrock ở dạng phân tích bằng:
đâu и , xác định ma trận .
Các phần tử khác XNUMX còn lại của ma trận bằng:
Ví dụ, trong không gian năm chiều N = 5, ma trận Hessian cho hàm Rosenbrock có dạng một dải:
Mã tính toán Hessian này cùng với mã để giảm thiểu hàm Rosenbrock bằng phương pháp gradient liên hợp (Newton):
def rosen_hess(x):
x = np.asarray(x)
H = np.diag(-400*x[:-1],1) - np.diag(400*x[:-1],-1)
diagonal = np.zeros_like(x)
diagonal[0] = 1200*x[0]**2-400*x[1]+2
diagonal[-1] = 200
diagonal[1:-1] = 202 + 1200*x[1:-1]**2 - 400*x[2:]
H = H + np.diag(diagonal)
return H
res = minimize(rosen, x0, method='Newton-CG',
jac=rosen_der, hess=rosen_hess,
options={'xtol': 1e-8, 'disp': True})
print(res.x)Optimization terminated successfully.
Current function value: 0.000000
Iterations: 24
Function evaluations: 33
Gradient evaluations: 56
Hessian evaluations: 24
[1. 1. 1. 0.99999999 0.99999999]Một ví dụ với định nghĩa về hàm tích Hessian và một vectơ tùy ý
Trong các vấn đề thực tế, việc tính toán và lưu trữ toàn bộ ma trận Hessian có thể yêu cầu tài nguyên bộ nhớ và thời gian đáng kể. Trong trường hợp này, thực sự không cần phải chỉ định chính ma trận Hessian, bởi vì quy trình tối thiểu hóa chỉ yêu cầu một vectơ bằng tích của Hessian với một vectơ tùy ý khác. Do đó, từ quan điểm tính toán, tốt hơn hết là xác định ngay một hàm trả về kết quả của tích Hessian với một vectơ tùy ý.
Hãy xem hàm hess, lấy vectơ cực tiểu hóa làm đối số đầu tiên và vectơ tùy ý làm đối số thứ hai (cùng với các đối số khác của hàm được giảm thiểu). Trong trường hợp của chúng tôi, việc tính tích Hessian của hàm Rosenbrock với một vectơ tùy ý không khó lắm. Nếu như p là một vectơ tùy ý thì tích giống như:
Hàm tính tích của Hessian và một vectơ tùy ý được truyền dưới dạng giá trị của đối số hessp cho hàm thu nhỏ:
def rosen_hess_p(x, p):
x = np.asarray(x)
Hp = np.zeros_like(x)
Hp[0] = (1200*x[0]**2 - 400*x[1] + 2)*p[0] - 400*x[0]*p[1]
Hp[1:-1] = -400*x[:-2]*p[:-2]+(202+1200*x[1:-1]**2-400*x[2:])*p[1:-1]
-400*x[1:-1]*p[2:]
Hp[-1] = -400*x[-2]*p[-2] + 200*p[-1]
return Hp
res = minimize(rosen, x0, method='Newton-CG',
jac=rosen_der, hessp=rosen_hess_p,
options={'xtol': 1e-8, 'disp': True})
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 24
Function evaluations: 33
Gradient evaluations: 56
Hessian evaluations: 66Thuật toán vùng tin cậy gradient liên hợp (Newton)
Điều kiện kém của ma trận Hessian và hướng tìm kiếm không chính xác có thể khiến thuật toán gradient liên hợp của Newton không hiệu quả. Trong những trường hợp như vậy, ưu tiên được dành cho (vùng tin cậy) liên hợp gradient Newton.
Ví dụ với định nghĩa ma trận Hessian:
res = minimize(rosen, x0, method='trust-ncg',
jac=rosen_der, hess=rosen_hess,
options={'gtol': 1e-8, 'disp': True})
print(res.x)Optimization terminated successfully.
Current function value: 0.000000
Iterations: 20
Function evaluations: 21
Gradient evaluations: 20
Hessian evaluations: 19
[1. 1. 1. 1. 1.]Ví dụ với hàm tích của Hessian và một vectơ tùy ý:
res = minimize(rosen, x0, method='trust-ncg',
jac=rosen_der, hessp=rosen_hess_p,
options={'gtol': 1e-8, 'disp': True})
print(res.x)Optimization terminated successfully.
Current function value: 0.000000
Iterations: 20
Function evaluations: 21
Gradient evaluations: 20
Hessian evaluations: 0
[1. 1. 1. 1. 1.]Phương pháp loại Krylov
Giống như phương pháp Trust-ncg, các phương pháp kiểu Krylov rất phù hợp để giải các bài toán quy mô lớn vì chúng chỉ sử dụng tích vectơ ma trận. Bản chất của chúng là giải quyết vấn đề trong vùng tin cậy bị giới hạn bởi không gian con Krylov bị cắt cụt. Đối với các bài toán không chắc chắn, tốt hơn nên sử dụng phương pháp này, vì nó sử dụng số lần lặp phi tuyến nhỏ hơn do số lượng tích vectơ ma trận trên mỗi bài toán con nhỏ hơn so với phương pháp Trust-ncg. Ngoài ra, lời giải của bài toán con bậc hai được tìm ra chính xác hơn so với việc sử dụng phương pháp Trust-ncg.
Ví dụ với định nghĩa ma trận Hessian:
res = minimize(rosen, x0, method='trust-krylov',
jac=rosen_der, hess=rosen_hess,
options={'gtol': 1e-8, 'disp': True})
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 19
Function evaluations: 20
Gradient evaluations: 20
Hessian evaluations: 18
print(res.x)
[1. 1. 1. 1. 1.]
Ví dụ với hàm tích của Hessian và một vectơ tùy ý:
res = minimize(rosen, x0, method='trust-krylov',
jac=rosen_der, hessp=rosen_hess_p,
options={'gtol': 1e-8, 'disp': True})
Optimization terminated successfully.
Current function value: 0.000000
Iterations: 19
Function evaluations: 20
Gradient evaluations: 20
Hessian evaluations: 0
print(res.x)
[1. 1. 1. 1. 1.]
Thuật toán giải gần đúng trong miền tin cậy
Tất cả các phương pháp (Newton-CG, Trust-ncg và Trust-krylov) đều rất phù hợp để giải quyết các vấn đề quy mô lớn (với hàng nghìn biến số). Điều này là do thực tế là thuật toán gradient liên hợp cơ bản bao hàm việc xác định gần đúng ma trận Hessian nghịch đảo. Giải pháp được tìm thấy lặp đi lặp lại mà không cần mở rộng rõ ràng Hessian. Vì bạn chỉ cần xác định hàm cho tích của một vectơ Hessian và một vectơ tùy ý, nên thuật toán này đặc biệt tốt khi làm việc với các ma trận thưa thớt (đường chéo dải). Điều này mang lại chi phí bộ nhớ thấp và tiết kiệm thời gian đáng kể.
Đối với các bài toán cỡ trung bình, chi phí lưu trữ và phân tích Hessian là không quan trọng. Điều này có nghĩa là có thể đạt được giải pháp với số lần lặp ít hơn, giải quyết gần như chính xác các bài toán con của vùng tin cậy. Để làm điều này, một số phương trình phi tuyến được giải lặp cho từng bài toán con bậc hai. Một giải pháp như vậy thường đòi hỏi 3 hoặc 4 phân tách Cholesky của ma trận Hessian. Kết quả là, phương pháp này hội tụ ít lần lặp hơn và yêu cầu tính toán hàm mục tiêu ít hơn so với các phương pháp vùng tin cậy được triển khai khác. Thuật toán này chỉ liên quan đến việc xác định ma trận Hessian hoàn chỉnh và không hỗ trợ khả năng sử dụng hàm tích của Hessian và một vectơ tùy ý.
Ví dụ về việc tối thiểu hóa hàm Rosenbrock:
res = minimize(rosen, x0, method='trust-exact',
jac=rosen_der, hess=rosen_hess,
options={'gtol': 1e-8, 'disp': True})
res.xOptimization terminated successfully.
Current function value: 0.000000
Iterations: 13
Function evaluations: 14
Gradient evaluations: 13
Hessian evaluations: 14
array([1., 1., 1., 1., 1.])Có lẽ chúng ta sẽ dừng ở đó. Trong bài viết tiếp theo, tôi sẽ cố gắng kể những điều thú vị nhất về cực tiểu hóa có điều kiện, ứng dụng cực tiểu hóa trong việc giải các bài toán gần đúng, cực tiểu hóa hàm một biến, cực tiểu hóa tùy ý và tìm nghiệm của hệ phương trình bằng cách sử dụng scipy.optimize bưu kiện.
Nguồn:
Nguồn: www.habr.com