Python에서 최소 제곱법을 호출하고 구현하는 방법

소위 선형 최소제곱법은 방정식 풀이의 연속으로 이해될 수 있습니다. 차이점은 미지의 양이 방정식의 수보다 훨씬 작을 때 풀 수 없는 문제가 발생한다는 것입니다. 최소제곱법의 핵심은 오차를 최소화하면서 알 수 없는 숫자에 값을 할당하는 것입니다.

최소제곱법은 매우 고전적인 알고리즘이며, 고등학교 때 이 이름을 접한 적이 있습니다. 나는 이전에 선형 최소 제곱의 원리에 대해 작성하고 이를 Python으로 구현했습니다: 최소 제곱 및 Python 구현, scipy에서 비선형 최소 제곱을 호출하는 방법: 비선형 최소 제곱 (문서 끝의 보충 내용) ; 및 희소 행렬에 대한 최소 제곱법: 희소 행렬 최소 제곱법.

다음은 numpy와 scipy에서 구현한 선형 최소제곱법을 설명하고, 둘의 속도를 비교합니다.

numpy 구현

최소 제곱법은 numpy에서 구현됩니다. 즉, a@x=b와 유사하게 x를 풀기 위해 lstsq(a,b)가 사용됩니다. 여기서 a는 M×N의 행렬이고, b는 다음과 같습니다. M개의 벡터 행은 선형 방정식 시스템을 푸는 것과 정확히 동일합니다. Ax=b와 같은 방정식 시스템의 경우 A가 전체 순위 시뮬레이션이면 x=A⁻¹b로 표현될 수 있고, 그렇지 않으면 x=(A^TA)^{로 표현될 수 있습니다. −1}A ^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 패키지

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

Supplement

Python은 비선형 최소 제곱법을 호출합니다

소개 및 생성자

scipy에서 비선형 최소 제곱법의 목적은 다음의 제곱합을 최소화하는 함수 집합을 찾는 것입니다. 오류 함수는 다음과 같은 공식으로 표현될 수 있습니다

여기서 ρ는 손실 함수를 나타내며, 이는 f_i(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为置信半径，假设在这个邻域内，目标函数可以近似为线性或二次函数，则可通过二次模型得到区间中的极小值点s_k。然后以这个极小值点为中心，继续优化信赖域所对应的区间。

雅可比矩阵

在了解了信赖域方法之后，就会明白雅可比矩阵在数值求解时的重要作用，而如何计算雅可比矩阵，则是接下来需要考虑的问题。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 중국어 웹사이트의 기타 관련 기사를 참조하세요!