希尔伯特矩阵的病态性问题

希尔伯特矩阵的病态性问题
（学号：200921100102 姓 名：陈丹丹 ）
一、问题叙述
希尔伯特矩阵是著名的病态矩阵。

n 阶Hilbert 矩阵为
⎥
⎥
⎥⎥⎦
⎤⎢
⎢⎢
⎢⎣⎡−++=)12/(1)1/(1/1)1/(13/12/1/12/11
n n n n n H n L L L L L L L 其条件数 Cond(H n ) 如下表所示
n
2 3 4 5 6 Cond(H n )
1.9281e+001 5.2406e+002 1.5514e+004 4.7661e+005 1.4951e+007
随着n 的增大，矩阵条件数迅速增加。

猜测：希尔伯特矩阵条件数以指数规律增长。

即，设矩阵阶数为n ，有
Cond(H n ) ≈ exp( a n + b )
对表中的条件数做对数变换，问题转化为线性拟合问题
ln[Cond(H n )] = a n + b ，( n = 2，3，4，5，6)
线性函数的两个系数分别为
a = 3.3935，
b = – 3.8811
故指数拟合函数为：
Cond(H n ) ≈ exp(3.3935 n –3.8811 )
拟合函数的残差向量为
r 2 r 3 r 4 r 5 r 6
9.9860e-001 -2.0231e+001 -6.8994e+002 -5.7825e+003 5.9012e+005 我们讨论以下三个问题：
（1）由于线性拟合所得残差向量较大，下面用实验确定：选择用几次多项式拟合较好
（2）设D 是由H n 的对角线元素开方构成的对角矩阵，令1
1
ˆ−−=D H D H n
n ，不难看出该矩阵仍然是对称正定矩阵，而且其对角元素全为1。

将H n 变换为n H ˆ的技术称为预处理（Preconditioning ）。

试分析函数)](/)ˆ(ln[n
n H Cond H Cond 随n 变化的规律，绘制该函数的曲线。

（3）如果只用左预处理，即D 取希尔伯特矩阵对角元构成对角矩阵，用D 的逆阵左乘H 后条件数如何。

二、问题分析
（1）从上面的表中我们可以看出用线性拟合得到的残差向量较大，从而我们就想分别用二到十次多项式做拟合，比较得出来的拟合函数的残差向量，看看几次多项式拟合的拟合函数的残差向量最小，从而几次多项式拟合较好。

（2）我们编写程序，画出函数)](/)ˆ(ln[n
n H Cond H Cond ，并通过函数曲线来观察它随n 的变化规律和函数值在零之下的n 的取值数。

(3)取D 为希尔伯特矩阵对角元构成的对角矩阵，我们考查用D 的逆阵左乘H 后条件数。

同时画))(/)((1
n n H cond H D cond In −的函数图，并与（2）中的函数图比较。

三、实验程序
MATLAB 程序段如下： 程序一：
C=[]; for k=2:6
H=hilb(k); %生成k 阶希尔伯特矩阵
C=[C,cond(H)]; %希尔伯特矩阵的条件数
end
LC=log(C);n=2:6; %对希尔伯特矩阵的条件数取对数 format short e
P=polyfit(n,LC,2) %二次多项式拟合 % P=polyfit(n,LC,3) %三次多项式拟合 % P=polyfit(n,LC,4) %四次多项式拟合 % P=polyfit(n,LC,5) %五次多项式拟合 % P=polyfit(n,LC,6) %六次多项式拟合 % P=polyfit(n,LC,7) %七次多项式拟合 % P=polyfit(n,LC,8) %八次多项式拟合 % P=polyfit(n,LC,9) %九次多项式拟合 % P=polyfit(n,LC,10) %十次多项式拟合 residure=C-exp(polyval(P,n)) %残差
程序二：
n=500; C0=[]; for k=2:n
H=hilb(k); %生成k 阶希尔伯特矩阵 D=diag(sqrt(diag(H)));
D1=inv(D); %生成预条件子
H1=D1*H*D1; %对希尔伯特矩阵进行预处理
C0=[C0,cond(H1)/cond(H)]; %预处理后的条件数与之前的条件数之比 end
LC0=log(C0); %取对数 k=2:n;
plot(k,LC0,'-') %画图 程序三：
n=500; C=[]; C1=[]; for k=2:n
H=hilb(k); %生成k 阶希尔伯特矩阵 D=diag(diag(H));
H1= inv(D) *H; %对希尔伯特矩阵进行左预处理 C=[C,cond(H)]; %计算原希尔伯特矩阵的条件数
C1=[C1,cond(H1)]; %计算左预处理后矩阵的条件数 end
LC=log(C); %取对数 LC1=log(C1); k=2:n;
plot(k,LC,'-',k,LC1, 'r:') %画图
程序四：
n=500; C0=[]; for k=2:n
H=hilb(k); %生成k 阶希尔伯特矩阵 D=diag(diag(H));
H1= inv(D)*H; %对希尔伯特矩阵进行左预处理
C0=[C0,cond(H1)/cond(H)]; %左预处理后的条件数与之前的条件数之比 end
LC0=log(C0); %取对数 k=2:n;
plot(k,LC0,'-') %画图
四、实验数据结果及分析
（1） 程序一的第八行分别用2:10次的多项式拟合，得到的拟合函数的残差向量如下表：
表1 拟合函数的残差向量
m
r 2 r 3 r 4
r 5
r 6
1 9.9860e-001 -2.0231e+001-6.8994e+002
-5.7825e+003
5.9012e+005 2 1.3285e-001 -7.7849e+000 4.2669e+001 5.2486e+003 -8.9912e+004 3 8.8489e-003 -9.6313e-001 4.2669e+001 -8.7592e+002
6.8615e+003 4 -1.0303e-013 0 -1.4552e-011 -1.6298e-009-
7.2643e-008 5 7.1054e-015 9.0949e-013 4.1837e-011 5.8208e-011 -3.3714e-007
6 -7.2831e-013 -7.5033e-012-4.1837e-011 -6.6939e-009
4.6007e-007
7 8.7041e-013 0
1.5098e-010 1.3621e-008 9.9093e-007
8 -7.7343e-012 9.3223e-012 -1.8008e-010 5.8208e-011 3.1684e-006 9 -2.7644e-011 1.8611e-010 9.2223e-010 3.0559e-008 -7.2438e-006
10 -5.3916e-010 -2.4954e-010
-4.1473e-009
-2.7642e-006
-2.0192e-004
其中i r 表示用多项式拟合))((i H cond In 的残差;m 表示用m 次多项式进行拟合。

我们从每行可以看出，随着阶数的增大，对希尔伯特矩阵条件数的对数进行多项式拟合的残差也处于增大的趋势。

再从每列，我们可以看出，在4=m 和5=m 的时候，残差是较小的。

从下面的图1更加直观地发现，用四次或者五次多项式对)6,,2())((L =i H cond In i 进行多项式拟合，得到的残差最小，也就是说拟合得精确度最好。

65432*r r r r r ：：：：：×+∇Ο
图1
（2） 在Matlab 中运行程序二，我们得到图2，如下：
图2
上图中的11ˆ1−−==D H D H H n
n n （H n 的对角线元素开方构成的对角矩阵），n 为希尔伯特
矩阵的阶数。

在图2中，给出了n 从2到500变化时，对应的函数)](/)ˆ(ln[n n H Cond H Cond 的变化曲线。

我们观察到，随希尔伯特矩阵阶数的增大，函数)](/)ˆ(ln[n n H Cond H Cond 在[-6,4]区间波动，主要集中在[-3,1]区间。

我们知道当)()ˆ(n n H Cond H Cond ≤时，有0)](/)ˆ(ln[≤n n H Cond H Cond ，在上图中，我们可以很容易的观察到，对于大部分n ，函数值都是小于或者等于零的，这说明n H 经过预处理后的n
H ˆ的条件数较小，在一定程度上改善了原希尔伯特矩阵的性质。

（3）下面给出)(n H cond 与)(1
n H D cond −（D 为希尔伯特矩阵对角元构成的对角矩阵）的对比图如下：
图3
从图3中，整体来说，对于大部分的n ，希尔伯特矩阵经过左预处理后的矩阵的条件数比原矩阵要小，但是在图中不是特别明显，下面我们给出))(/)((1
n n H cond H D cond In −的函数图，如下：
图4
从图4中，可以观察到函数值在区间[-6，6]波动，主要集中在区间[-3,2]，对于部分的n，进行左预处理后得到的矩阵的条件数降低了，然而对于另一部分的n，则条件数增大了。

将它与图2进行对比，我们发现，图2中函数值小于等于零所对应的n的取值个数要多些，这也就说明了，对原希尔伯特矩阵进行左右预处理技术对于改善原希尔伯特矩阵的性质要优于对希尔伯特矩阵进行左预处理的技术。

五、实验结论
(1) 从表1和图1中，有m取4和5时，拟合函数的残差向量较小，也就是说，四次多项式
拟合和五次多项式拟合较好。

(2) 由于希尔伯特矩阵病态的性质，对于对原希尔伯特矩阵进行预处理后的新希尔伯特矩阵
的条件数在一定的区间呈波动的变化规律。

从整体上观察，对于大多数的n，进行预处理后的希尔伯特矩阵的条件数有明显的降低，这就说明预处理改善了大多数阶数的希尔伯特矩阵的性质。

(3) 对原希尔伯特矩阵进行左预处理虽然对于部分n，使左预处理后的矩阵的条件数较原矩
阵有所降低，但是还有很大部分的n所对应的希尔伯特矩阵，经过这样的左预处理后，反而得到了更大条件数的矩阵，这跟希尔伯特本身的病态性有着紧密的联系。

我们还得到，左预处理技术作用于希尔伯特矩阵改善原希尔伯特矩阵性质的效果不如左右预处理技术。