소위 선형 최소제곱법은 방정식 풀이의 연속으로 이해될 수 있습니다. 차이점은 미지의 양이 방정식의 수보다 훨씬 작을 때 풀 수 없는 문제가 발생한다는 것입니다. 최소제곱법의 핵심은 오차를 최소화하면서 알 수 없는 숫자에 값을 할당하는 것입니다.
최소제곱법은 매우 고전적인 알고리즘이며, 고등학교 때 이 이름을 접한 적이 있습니다. 나는 이전에 선형 최소 제곱의 원리에 대해 작성하고 이를 Python으로 구현했습니다: 최소 제곱 및 Python 구현, scipy에서 비선형 최소 제곱을 호출하는 방법: 비선형 최소 제곱 (문서 끝의 보충 내용) ; 및 희소 행렬에 대한 최소 제곱법: 희소 행렬 최소 제곱법.
다음은 numpy와 scipy에서 구현한 선형 최소제곱법을 설명하고, 둘의 속도를 비교합니다.
최소 제곱법은 numpy에서 구현됩니다. 즉, a@x=b와 유사하게 x를 풀기 위해 lstsq(a,b)가 사용됩니다. 여기서 a는 M×N의 행렬이고, b는 다음과 같습니다. M개의 벡터 행은 선형 방정식 시스템을 푸는 것과 정확히 동일합니다. Ax=b와 같은 방정식 시스템의 경우 A가 전체 순위 시뮬레이션이면 x=A−1b로 표현될 수 있고, 그렇지 않으면 x=(ATA)로 표현될 수 있습니다. −1A Tb.
b가 M×K의 행렬이면 각 열에 대해 x 집합이 계산됩니다.
4개의 반환 값이 있는데, 피팅으로 얻은 x, 피팅 오류, 행렬 a의 순위, 행렬 a의 단일 값 형식입니다.
import numpy as np np.random.seed(42) M = np.random.rand(4,4) x = np.arange(4) y = M@x xhat = np.linalg.lstsq(M,y) print(xhat[0]) #[0. 1. 2. 3.]
scipy.linalg도 함수 이름은 lstsq이고 해당 매개변수 목록은
lstsq(a, b, cond=None, overwrite_a=False, overwrite_b=False, check_finite=True, lapack_driver=None)
입니다. 여기서 a와 b는 모두 재정의 가능한 스위치를 제공합니다. 또한 이 함수는 linalg의 많은 함수에 있는 옵션인 유한성 검사도 지원합니다. 반환 값은 numpy의 최소 제곱 함수와 동일합니다.
cond는 특이값 임계값을 나타내는 부동 소수점 매개변수입니다. 특이값이 cond보다 작으면 폐기됩니다.
lapack_driver는 LAPACK에서 어떤 알고리즘 엔진이 선택되었는지 나타내는 문자열 옵션이며 선택적으로 'gelsd', 'gelsy', 'gelss'입니다.
import scipy.linalg as sl xhat1 = sl.lstsq(M, y) print(xhat1[0]) # [0. 1. 2. 3.]
마지막으로 두 세트의 최소 제곱 함수의 속도를 비교합니다
from timeit import timeit N = 100 A = np.random.rand(N,N) b = np.arange(N) timeit(lambda:np.linalg.lstsq(A, b), number=10) # 0.015487500000745058 timeit(lambda:sl.lstsq(A, b), number=10) # 0.011151800004881807
이번에는 행렬 차원을 500으로 확대해도 둘 사이에 큰 차이가 없습니다. 2개도 대여섯 개예요.
N = 500 A = np.random.rand(N,N) b = np.arange(N) timeit(lambda:np.linalg.lstsq(A, b), number=10) 0.389679799991427 timeit(lambda:sl.lstsq(A, b), number=10) 0.35642060000100173
Python은 비선형 최소 제곱법을 호출합니다
소개 및 생성자
scipy에서 비선형 최소 제곱법의 목적은 다음의 제곱합을 최소화하는 함수 집합을 찾는 것입니다. 오류 함수는 다음과 같은 공식으로 표현될 수 있습니다
여기서 ρ는 손실 함수를 나타내며, 이는 fi(x)의 전처리로 이해될 수 있습니다.
scipy.optimize는
least_squares(fun, x0, jac, bounds, method, ftol, xtol, gtol, x_scale, f_scale, loss, jac_sparsity, max_nfev, verbose, args, kwargs)
로 정의되는 비선형 최소제곱 함수인 less_squares를 캡슐화합니다. 그 중 func와 x0은 필수 매개변수이고, func는 풀려는 함수이고, x0은 함수 입력의 초기값입니다. is no 기본값은 반드시 입력해야 하는 매개변수입니다.
bound는 솔루션 간격이며 기본값은 (−∞,∞)입니다. verbose가 1이면 종료 출력이 발생하고 verbose가 2이면 작업 중 추가 정보가 인쇄됩니다. 또한 다음과 같은 매개변수를 사용하여 오류를 제어하는데 비교적 간단합니다.
| 기본값 | Remarks | |
|---|---|---|
| ftol | 10-8 | 기능 공차 |
| xtol | 1 0-8 | 독립변수 공차 |
| gtol | 10-8 | 경사 허용오차 |
| x_scale | 1.0 | 변수의 특성 척도 |
| f_scale | 1.0 | Mar 잔차의 최종값 |
loss为损失函数，就是上面公式中的ρ \rhoρ，默认为linear，可选值包括
迭代策略
上面的公式仅给出了算法的目的，但并未暴露其细节。关于如何找到最小值，则需要确定搜索最小值的方法，method为最小值搜索的方案，共有三种选项，默认为trf
trf：即Trust Region Reflective，信赖域反射算法
dogbox：信赖域狗腿算法
lm：Levenberg-Marquardt算法
这三种方法都是信赖域方法的延申，信赖域的优化思想其实就是从单点的迭代变成了区间的迭代，由于本文的目的是介绍scipy中所封装好的非线性最小二乘函数，故而仅对其原理做简略的介绍。
其中r为置信半径，假设在这个邻域内，目标函数可以近似为线性或二次函数，则可通过二次模型得到区间中的极小值点sk。然后以这个极小值点为中心，继续优化信赖域所对应的区间。
雅可比矩阵
在了解了信赖域方法之后，就会明白雅可比矩阵在数值求解时的重要作用，而如何计算雅可比矩阵，则是接下来需要考虑的问题。jac参数为计算雅可比矩阵的方法，主要提供了三种方案，分别是基于两点的2-point；基于三点的3-point；以及基于复数步长的cs。一般来说，三点的精度高于两点，但速度也慢一倍。
此外，可以输入自定义函数来计算雅可比矩阵。
测试
最后，测试一下非线性最小二乘法
import numpy as np
from scipy.optimize import least_squares
def test(xs):
_sum = 0.0
for i in range(len(xs)):
_sum = _sum + (1-np.cos((xs[i]*i)/5)*(i+1))
return _sum
x0 = np.random.rand(5)
ret = least_squares(test, x0)
msg = f"最小值" + ", ".join([f"{x:.4f}" for x in ret.x])
msg += f"\nf(x)={ret.fun[0]:.4f}"
print(msg)
'''
最小值0.9557, 0.5371, 1.5714, 1.6931, 5.2294
f(x)=0.0000
'''위 내용은 Python에서 최소 제곱법을 호출하고 구현하는 방법의 상세 내용입니다. 자세한 내용은 PHP 중국어 웹사이트의 기타 관련 기사를 참조하세요!
![[웹 프런트엔드] Node.js 빠른 시작](https://img.php.cn/upload/course/000/000/067/662b5d34ba7c0227.png)
