x 1 、 x 2 、 x 3 x_{1}、x_{2}、x_{3} x1、x2、x3为待定参数,为了确定这三个参数,测得 t 1 、 t 2 、 y t_{1}、t_{2}、y t1、t2、y的5组数据:
y = x 1 x 3 t 1 1 + x 1 t 1 + x 2 t 2 y= \frac{ x_{1}x_{3}t_{1} } { 1+x_{1}t_{1}+x_{2}t_{2} } y=1+x1t1+x2t2x1x3t1
t 1 t_{1} t1 | t 2 t_{2} t2 | y y y |
---|---|---|
1 | 1 | 0.126 |
2 | 1 | 0.219 |
1 | 2 | 0.076 |
2 | 2 | 0.126 |
0.1 | 0 | 0.186 |
原问题的最小二乘模型
m
i
n
f
(
x
)
T
f
(
x
)
minf(x)^Tf(x)
minf(x)Tf(x)
f
(
x
)
=
[
f
1
(
x
)
f
2
(
x
)
f
3
(
x
)
f
4
(
x
)
f
5
(
x
)
]
=
[
(
x
1
x
3
1
+
x
1
+
x
2
−
0.126
)
2
(
2
x
1
x
3
1
+
2
x
1
+
x
2
−
0.219
)
2
(
x
1
x
3
1
+
x
1
+
2
x
2
−
0.076
)
2
(
2
x
1
x
3
1
+
x
1
+
2
x
2
−
0.126
)
2
(
0.1
x
1
x
3
1
+
x
1
−
0.186
)
2
]
f(x)=\left[\begin{array}{l} f_{1}(x) \\ f_{2}(x) \\ f_{3}(x) \\ f_{4}(x) \\ f_{5}(x) \end{array}\right]=\left[\begin{array}{l} \left(\frac{x_{1} x_{3}}{1+x_{1}+x_{2}}-0.126\right)^{2} \\ \left(\frac{2 x_{1} x_{3}}{1+2 x_{1}+x_{2}}-0.219\right)^{2} \\ \left(\frac{x_{1} x_{3}}{1+x_{1}+2 x_{2}}-0.076\right)^{2} \\ \left(\frac{2 x_{1} x_{3}}{1+x_{1}+2 x_{2}}-0.126\right)^{2} \\ \left(\frac{0.1 x_{1} x_{3}}{1+x_{1}}-0.186\right)^{2} \end{array}\right]
f(x)=⎣⎢⎢⎢⎢⎡f1(x)f2(x)f3(x)f4(x)f5(x)⎦⎥⎥⎥⎥⎤=⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎡(1+x1+x2x1x3−0.126)2(1+2x1+x22x1x3−0.219)2(1+x1+2x2x1x3−0.076)2(1+x1+2x22x1x3−0.126)2(1+x10.1x1x3−0.186)2⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎤
那么如果选用高斯方法,则迭代公式为
x
k
+
1
=
x
k
−
(
A
k
T
A
k
)
−
1
A
k
T
f
k
x_{k+1}=x_{k}-\left(A_{k}^{T} A_{k}\right)^{-1} A_{k}^{T} f_{k}
xk+1=xk−(AkTAk)−1AkTfk
其中,
A
(
x
k
)
=
[
∇
f
1
(
x
k
)
T
⋮
∇
f
m
(
x
k
)
T
]
=
[
∂
f
1
(
x
k
)
∂
x
1
⋯
∂
f
1
(
x
k
)
∂
x
n
⋯
⋯
⋯
∂
f
m
(
x
k
)
∂
x
1
⋯
∂
f
m
(
x
k
)
∂
x
n
]
A\left(x_{k}\right)=\left[\begin{array}{c} \nabla f_{1}\left(x_{k}\right)^{T} \\ \vdots \\ \nabla f_{m}\left(x_{k}\right)^{T} \end{array}\right]=\left[\begin{array}{ccc} \frac{\partial f_{1}\left(x_{k}\right)}{\partial x_{1}} & \cdots & \frac{\partial f_{1}\left(x_{k}\right)}{\partial x_{n}} \\ \cdots & \cdots & \cdots \\ \frac{\partial f_{m}\left(x_{k}\right)}{\partial x_{1}} & \cdots & \frac{\partial f_{m}\left(x_{k}\right)}{\partial x_{n}} \end{array}\right]
A(xk)=⎣⎢⎡∇f1(xk)T⋮∇fm(xk)T⎦⎥⎤=⎣⎢⎡∂x1∂f1(xk)⋯∂x1∂fm(xk)⋯⋯⋯∂xn∂f1(xk)⋯∂xn∂fm(xk)⎦⎥⎤
等价问题的最小二乘模型
对于线性问题
A
X
=
b
AX=b
AX=b,最小二乘解为
X
=
(
A
T
A
)
−
1
A
T
b
X=(A^{T}A)^{-1}A^{T}b
X=(ATA)−1ATb.
由于是非线性方程,可以通过一定的转换将其变为线性模型,具体操作如下,将该算式两边同时倒数,则有:
1 y = 1 t 1 ∗ 1 x 1 x 3 + t 2 t 1 ∗ x 2 x 1 x 3 + 1 x 3 \frac{1}{y}= \frac {1}{t_{1}}*\frac {1}{x_{1}x_{3}}+ \frac {t_2}{t_{1}}*\frac {x_{2}}{x_{1}x_{3}}+ \frac{1}{x_{3}} y1=t11∗x1x31+t1t2∗x1x3x2+x31
如果令
Y
=
1
y
,
X
1
=
1
x
1
x
3
,
X
2
=
x
2
x
1
x
3
,
X
3
=
1
x
3
Y=\frac{1}{y},X_{1}=\frac {1}{x_{1}x_{3}},X_{2}=\frac {x_{2}}{x_{1}x_{3}},X_{3}=\frac{1}{x_{3}}
Y=y1,X1=x1x31,X2=x1x3x2,X3=x31
A
1
=
1
t
1
,
A
2
=
t
2
t
1
A_{1}=\frac {1}{t_{1}},A_{2}=\frac{t_{2}}{t_{1}}
A1=t11,A2=t1t2
那么原方程可以转换为:
(
A
1
,
A
2
,
1
)
∗
(
X
1
,
X
2
,
X
3
)
T
=
Y
(A_{1},A_{2},1)*(X_{1},X_{2},X_{3})^{T}=Y
(A1,A2,1)∗(X1,X2,X3)T=Y
若把
Y
Y
Y看成
b
b
b,则有:
A
X
=
b
AX=b
AX=b
因此,其最小二乘法的数学模型为:
m
i
n
(
A
1
X
1
+
A
2
X
2
+
X
3
−
Y
)
2
min(A_{1}X_{1}+A_{2}X_{2}+X_{3}-Y)^{2}
min(A1X1+A2X2+X3−Y)2
A
A
A矩阵的维数是5×3
,
A
T
A
A^TA
ATA 的维数是3×3
,
A
T
b
A^Tb
ATb的维数是3×1
那么代入
X
=
(
A
T
A
)
−
1
A
T
b
X=(A^{T}A)^{-1}A^{T}b
X=(ATA)−1ATb可以求出来
X
1
、
X
2
、
X
3
X_{1}、X_{2}、X_{3}
X1、X2、X3.
很容易推导出
x
1
=
X
3
X
1
、
x
2
=
X
2
X
1
、
x
3
=
1
X
3
x_{1}=\frac {X_{3}}{X_{1}}、x_{2}=\frac{X_{2}}{X_{1}}、x_{3}=\frac{1}{X_{3}}
x1=X1X3、x2=X1X2、x3=X31
编写matlab程序如下:
A=[1 1 1;0.5 0.5 1;1 2 1;0.5 1 1;10 0 1;]
y=[0.126,0.219,0.076,0.126,0.186]
t1=[1,2,1,2,0.1]
t2=[1,1,2,2,0]
Y=1./y
b=Y'
X=inv(A'*A)*A'*b
x1=X(3)/X(1)
x2=X(2)/X(1)
x3=1/X(3)
for i = 1:5
y_fit(i)=(t1(i)*x1*x3)/(t1(i)*x1 + t2(i)*x2 + 1)
end
error = sum (y_fit - y)^2
结果如下:
y_fit = [0.1283 0.2054 0.0751 0.1312 0.1860]
error=4.9186e-05