Transcript E:WHU教学2010Lecture2008Lec10_LSPLec10_LSP

Least Squares Problems
问题的提出
插值法是使用插值多项式来逼近未知或复杂函数的，它 要求插值函
数与被插函数在插值节点上函数值相同 ，而在其他点上没有要求。在非
插值节点上有时函数值会相差很大 。若要求在被插函数的定义区间上，
所选近似函数都能与被插函数有较好的近似，就是最佳逼近问题。
最佳逼近是在函数空间 M中选 P(x) 满足
max
a xb
f ( x )  p( x )  min 
(*)
但由于绝对值函数不宜进行分析运算,常将上式化为

b
a
 ( x) (
f ( x)  p( x))
2
dx  min
来讨论 ，于是最佳逼近问题变为最佳平方逼近问题 ,而离散的最佳平方逼
进问题就是常说的曲线拟合
  ( f ( xi )  p( xi ))
m
i 0
2
i
它们都可用最小二乘法求解。
 min
曲线拟合的最小二乘法
最小二乘原理
 ( x)
当由实验提供了大量数据时,不能要求拟合函数
在数据点( xi , yi )
处的偏差,即
 i   ( xi )  yi
(i=1,2,…,m) 严格为零,但为了使近似曲线尽量反映所给数
据点的变化趋势 ,需对偏差有所要求.通常要求偏差平方和
m
|
i 1
2
i
|   ( ( x i )  y
m
i 1
最小,此即称为最小二乘原理
2
i)
最小二乘法的求法
设近似方程为 : (共有m组数据且m  n)
y* 
a
0
0

( x )  a1
1
( x )  ...  a n
m
 (a0 , a1 ,..., a n )   ( y i 
*
i 1
y)
2
i

n
( x)
 min
对函数求偏导数并令其为零, 可得 :

0
aj
 2(  a  ( x )  y )  ( x )  0
m
n
i 1
k 0
k
 m n
得 : 2   a k
 i 1 k  0
若引入记号 :
i
k



j
j
k
,
( xi )



i
j
i
m
j
( xi )  
i 1
y
   ( x )
i

( xi )   0
j

m
k
, f 
i
j
i 1
 ( x ) y
m
i 1
j
i
i
k
( xi )
n
 ,
则有 :  ak 

k 0
j

  
k

, f 
j

( j  0,1,..., n)
 可得矩阵







 , 
 , 
 , 
 , 
...
...
...
...
0
0
1
0
0
1
1
1
 ,   , 
n
0
n
1
...
...
 ,  a    , f 
 ,   a     , f 
0
n
0
1
  ... 
...
 
,
 a 
n
n   n
1
n
  
0






,f 
n

1

...

可知
当

0
( x ),

1
( x ),...

n
( x ) 线性无关时 存在唯一解
,
ai ( i  0,1,..., n)
n
a 
i 0
i
i
( x )就是所求的拟合函数
用矩阵表达：
给定一组测量或观测值
yi  y(ti ), i  1,..., m
假设y是n个函数的线性组合
y(t )  11 (t )     nn (t )
设m×n矩阵X定义为
xij   j (ti )
所以
y  Xβ
min X β  y
R( X)   N ( XT )
β
XT X β  XT y
1、求逆的过程困难
T
2
2、 cond ( X X )  (cond ( X ))
1 1


X   0 
0  


2

1


1 
T

X X  
2
1  
 1
作为曲线拟合的一种常见情况, 拟合函数常为代数多项式
:,
即拟合函数 :
 ( x)  a0  a1 x  a2 x  ...  an x
2
n
 0  1 x1   2 x12     m x1m  y1

2
m
 0  1 x2   2 x2     m x2  y2


   x   x 2     x m  y
1 n
2 n
m n
n
 0
写成矩阵形式为：
其中，
1

1

X


1
Xβ  y
 0 
x12  x1m 
 
2
m 
x 2 x 2  x2 
1 

, β
,

 
   

 
xn xn2  xnm 
 m 
x1
 y1 
y 
2

y
 
 
 yn 
XT Xβ  y

 n

 n
xi
T


X X  i 1

 n

m
  xi
 i 1
n
x
i
i 1
n
x
2
i
i 1

n
x
i 1
m 1
i
n



m
  yi 
  xi 
 i 1

i 1

n
n



m 1
  xi  XT y    xi yi 
,
i 1
 i 1


 
 n  
n
 xm y 
2m 
  xi 
 i i 
 i 1

i 1

n
y ( x)   0  1 x

 n
XT X   n

  xi
 i 1

xi 

i 1

n
2
xi 

i 1

n
 n

  yi 

XT y   ni 1


  xi yi 
 i 1

9
8
7
6
5
4
3
2
1
1
2
3
4
5
6
7
8
9
10
直线：
y(t )  1t   2
多项式：  j (t )  t n j , j  1,..., n
y (t )  1t n 1  ...   n 1t   n
t n j
有理函数：  j (t ) 
1t n 1  ...   n 1t   n
1t n 1  ...   n 1t  
y (t ) 
1t n 1  ...   n 1t   n
指数：
 j (t )  e
 jt
y (t )  1e  1t  ...   n e  nt
Log-linear:
y (t )  Ke
t
log y(t )  1t   2
Gaussians:
 j (t )  e
 t  j

 j

y (t )  1e




2
 t  1 


 1 
2
 ...   n e
 t n

 n



2
评注:最小二乘法方曲线拟和是实验数据处理的常用方法。最佳平方逼近可以在一个区间上
比较均匀的逼近函数且具有方法简单易行，实效性大，应用广泛等特点。但当正规方
程阶数较高时，往往出现病态。因此必须谨慎对待和加以巧妙处理。有效方法之一是
引入正交多项式以改善其病态性。
正交多项式
付立叶级数，函数系
1,cosx,sinx,cos2x,sin2x,…cosnx,sinnx,…
中，由于任意两个函数乘积在区间[-，+]上的积分都
等于零，则说这个函数系在[-，+]上是正交的，并称
这个函数系为正交函数系。
下面给出正交函数系
定义：设函数f(x),g(x)[a,b],且
( f , g) 
  ( x ) f ( x ) g ( x )dx  0
b
a
则称f(x)与g(x)在[a,b]上带权(x)正交,
在[a,b]上连续的函数0(x), 1(x), 2(x),... k(x)...,
满足
b
jk
0
( j ,  k )   ( x ) j ( x ) k ( x )d ( x )  
a
jk
 Ak  0

则称该函数系是在区间[a,b]上带权(x)正交函数系.
下面介绍与上述定义有关的几个概念，然后引出正交多项
的概念，最后再介绍正交多项式的性质以及几种常见的正
交多项式。
1.权函数:(1)设[a,b]是有限或无限区间, (x)是定义在[a,b]
上的非零可积函数，若其满足
b

(1)  ( x )dx  0
a
b

(2) x n  ( x )dx
存在 n  1,2,
a
则称(x)是[a,b]上的一个权函数。
2 内积与范数
设f(x),g(x)[a,b], (x)是[a,b]上的一个权函数，称
( f , g) 
  ( x ) f ( x ) g ( x )dx
b
a
为f(x)与g(x)在为 [a,b]上以权函数(x)的内积。
显然，对于任意实数a,b,有
( f , ag  bh)  a( f , g)  b( f , h)
称
f
1
2
2

b
 ( f , f )  (  ( x) f 2 ( x)dx)
a
为f(x)的带权(x)的2—范数。
1
2
3． 正交多项式。最高幂项的系数a n  0 的几次多项式 g n ( x ), n  0,1, …若满足
0, i  j

( gi , g j )  
则g n (x ) 称为在[a,b]上带权  (x ) 正交，g n (x ) 称为[a,b]上
2
 g j  0, i  j
带权  (x ) 的几次正交多项式。
q( x), g n ( x))  0
正交多项式的性质
定理1 [a,b]上带权(x)的正交多项式系{gn(x)}一定是
[a,b]
上线性无关的函数系。
定理2 设是{gn(x)}[a,b]上带权(x)的正交多项式系，则对于任
何次数不高于n-1的多项式q(x),总有
(q(x), gn(x))=0
( n=1,2,…）
定理3 n次正交多项式gn(x)有n个互异定根，且全部若在(a,b)内。
定理4:任何相邻的三个正交多项式,都具有下列递推关系式
gn+1(x)=(nx-n)gn(x)-n-1gn-1(x)
常见的正交多项式
•勒让德多项式(Legendre)
•切比雪夫多项式(Chebyshev)
•拉盖尔多项式(Laguerre)
•埃尔米特多项式 (Hermite)
勒让德多项式(Legendre)
[-1,1] , (x)=1
递推关系:
P0(x)=1, P1(x)=x,
P3 ( x ) 
Pn 1 ( x) 
1
2
2 n 1
n 1
1 dn
Pn ( x )  n  n ( x 2  1) 2
2 ! dx
P2 ( x )  12 (3x 2  1)
(5 x 3  3 x )
xPn ( x)  nn1 Pn 1 ( x), n  1,2,3..... n
切比雪夫多项式(Chebyshev)
[1,1],  ( x ) 
1
1 x
2
,
Tn(x)=cos(narccosx)
递推关系:
T0(x)=1 , T1(x)=x , T2(x)=2x2-1 ,
T3(x)= 4x3-3x,………
T0 ( x )  1, T1 ( x )  x


Tn 1 ( x )  2 x, Tn ( x )  Tn 1 ( x ), n  1,2,3.... n
拉盖尔多项式(Laguerre)
[0,+), (x)=e-x
n
d
Ln ( x )  e x  n ( x n  e  x )
dx
埃尔米特多项式 (Hermite)
(- ,+), (x)=e-x2
n
d
(n)
 x2
H n ( x )  ( 1) e  n e
dx
x2
举例：
中国人口预测
y  1t 3   2t 2   3t   4
s  (t  1979) / 30
y  1s 3   2 s 2   3 s   4
p = 1e4*[ ]';
t = (1949:1:2005)';
x = (1940: 0.2 :2019)';
w = 2010;
d = 3;
s = (t-1979)/30; c = polyfit(s,p,d);
s = (x-1979)/30; y = polyval(c,s);
s = (w-1979)/30; z = polyval(c,s);
plot(t,p,'o',x,y,'b-',w,z,'r*')
s = (t-1979)/30; c = polyfit(s,p,3);
X = [s.^3, s.^2, s, ones(length(s),1)]; X\p
Matlab中避免用法方程求解
X\y
QR分解：
X  QR
Q正交阵(orthogonal)
QT Q  I
R上三角阵
求解： R  QT y
Householder 变换
Xβy
H n  H 2 H1X  R
Rβz
H n H 2 H1y  z
   || x ||,
u  x   ek ,
  2 / || u ||2  1 /  u k ,
H  I  uu T .
1950
55196
1960
66207
1970
82992
1980
98705
1990
114333
2000
126743
s:
-0.9667
-0.6333
-0.3000
0.0333
0.3667
0.7000
y:
X:
0.9344 -0.9667 1.0000
0.4011 -0.6333 1.0000
0.0900 -0.3000 1.0000
0.0011 0.0333 1.0000
0.1344 0.3667 1.0000
0.4900 0.7000 1.0000
s  (t  1979) / 30
y  1s 2   2 s1   3
55196
66207
82992
98705
114333
126743
s = ( (1950:10:2000)'-1979 )/30;
X = [s.^2, s, ones(size(s),1)]
y=[]
[m,n] = size(X);
for k = 1:min(m-1,n)
i = k:m;
u = X(i,k);
sigma = norm(u);
if sigma ~= 0
if u(1) ~= 0, sigma = sign(u(1))*sigma;
end
u(1) = u(1) + sigma;
rho = 1/( sigma *u(1));
u(1) = u(1) + sigma;
rho = 1/( sigma *u(1));
X(i,k) = 0;
X(k,k) = -sigma;
j = k+1:n;
v = rho*(u'* X(i,j));
X(i,j) = X(i,j) - u*v;
if ~isempty(y)
tau = rho*(u'* y(i));
y(i) = y(i) - tau*u;
end
end
end
beta = X(1:n,1:n)\y(1:n)
polyval(beta,(2010-1979)/30)
X=
-1.1403 0.6945
0
-0.3122
0
-0.2279
0
0.0342
0
0.4743
0
1.0923
X=
-1.1403
0
0
0
0
0
X=
-1.1403
0
0
0
0
0
-1.7987
0.4589
0.8786
0.9985
0.8186
0.3390
y = 1.0e+005 *
-1.4311
0.2787
0.7439
0.9860
1.0148
0.7991
0.6945 -1.7987
1.2525 0.3587
0
0.8640
0
1.0007
0
0.8490
0
0.4090
y = 1.0e+005 *
-1.4311
0.9033
0.8349
0.9723
0.8255
0.3630
0.6945 -1.7987
1.2525 0.3587
0
-1.6236
0
0
0
0
0
0
y = 1.0e+005 *
-1.4311
0.9033
-1.5667
0.0062
0.0058
-0.0318
β  (XT X) 1 XT y
XT Xβ  XT y
AXA  A
广义逆
其他广义逆: Group inverse,
Drazin inverse,…
XAX  X
AX 
XA T
T
 AX
 XA

列满秩
A   AT A
是左逆
A A  I
不是右逆,但满足
C
A  P
0

1
X
1 / x, x  0
x 
x0
 0,

C 0 1
AD  P 
P

 0 0
AT
min AX  I
0  1
P

N
F
秩亏损(Rank Deficiency)
接近秩亏损
X=
1
4
7
10
13
y=
2 3
5 6
8 9
11 12
14 15
16
17
18
19
20
X\y
pinv(X)*y
X\y-pinv(X)*y
ans =
ans =
ans =
-7.5000
0
7.8333
-7.5556
0.1111
7.7778
0.0556
-0.1111
0.0556
非线性最小二乘问题
y  1e1t   2e 2t
lsline
lsqlin.
lsqnonneg
lsqnonlin,
lsqcurvefit
fminunc
fmincon
fminsearch