我第一次遇到希尔伯特矩阵当我在做个人研究下教授约翰·托德加州理工学院在1960年。它一直是我职业生涯的一部分。

内容

大卫希尔伯特

在20世纪,大卫希尔伯特是世界上最著名的数学家。他介绍了矩阵,现在以他的名字命名在1895年的一篇论文中。矩阵的元素,这是连续正整数的倒数,沿着antidiagonals是不变的。

$ $ h_ {i, j} = \压裂{1}{我+ j - 1}, \ \我,j = 1: n $ $

格式老鼠H5 = hilb (5)

H5 = 1 1/2 1/3 1/4 1/2 1/3 1/4 1/5 1/5 1/6 1/3 1/4 1/5 1/6 1/7 1/5 1/6 1/7 1/8 1/5 1/6 1/7 1/8 1/9 1/4

乳胶(信谊(H5));

$ $ H_5 = \离开(\开始{数组}{ccccc} 1 & \压裂{1}{2}& \压裂{1}{3}& \压裂{1}{4}& \压裂{1}{5}\ \ \压裂{1}{2}& \压裂{1}{3}& \压裂{1}{4}& \压裂{1}{5}& \压裂{1}{6}\ \ \压裂{1}{3}& \压裂{1}{4}& \压裂{1}{5}& \压裂{1}{6}& \压裂{1}{7}\ \ \压裂{1}{4}& \压裂{1}{5}& \压裂{1}{6}& \压裂{1}{7}& \压裂{1}{8}\ \ \压裂{1}{5}& \压裂{1}{6}& \压裂{1}{7}& \压裂{1}{8}& \压裂{1}{9}\结束数组{}\右)$ $

这是一幅画。连续生成表面非常光滑。

H12 = hilb (12);冲浪(日志(H12))视图(60 60)

最小二乘

希尔伯特这个矩阵很感兴趣,因为它出现在一个连续函数的最小二乘逼近多项式的单位间隔,表达了在传统基础上的单项的线性组合。

$ $ p (x) = \ sum_ {j = 1} ^ n c_j x ^ {j - 1} $ $

正常的方程的系数矩阵元素

$ $ \ int_0 ^ 1 x ^ {i + j2} dx \ = \ \压裂{1}{i + j - 1} $ $

属性

一个希尔伯特矩阵有许多有用的特性。

希尔伯特矩阵的每一列几乎是其他列的倍数。所以几乎是线性相关的列和矩阵是接近,但不完全是,奇异。

hilb

MATLAB总是有功能hilb和invhilb计算希尔伯特矩阵及其逆。的身体hilb目前仅两行。

J = 1: n;H = 1. / (J ' + J - 1);

我们经常提到这是一个很好的例子单例扩张。一个列向量被添加到一个行向量产生一个矩阵,然后减去一个标量的矩阵,最后元素的倒数产生结果。

逆

可以表达希尔伯特的逆矩阵元素的二项式系数。原因我已经忘记,我总是使用新台币对美元H ^ {1} $。

$ $ t_ {i, j} =识别(1)^ {i + j} (i + j - 1) {{n +张}\选择{n-j}} {{n + j - 1} \选择{n}} {{i + j2} \选择张{}}$ $ ^ 2

希尔伯特逆矩阵的元素是整数,但它们大整数。最大的对角线上。为订单13日最大的元素

$ $ (T_{13}识别)_ {9,9}\ = \ 100863567447142500 $ $

这是在美元10 ^{17}$和大于双精度flintmax。

格式再有flintmax = flintmax

flintmax = 9.007199254740992 e + 15

因此,它可以代表最大的元素只有在n \ le 12美元。

帮助条目invhilb包括一个句子继承了我最初的Fortran MATLAB。

N小于15的结果是准确的。

这是误导。它应该说

N < = 12的结果是准确的。

我申请一个错误报告。

invhilv

的身体invhilb首先设置p对订单n。双嵌套为递归循环然后评估二项式系数公式,避免不必要的整数溢出。

dbtypeinvhilb18:28

18 p = n;19因为我= 1:n 20 r = p * p;21 H(我)= r /(2张);22 j = i + 1: n 23 r = - ((n-j + 1) * r * (n + j - 1)) / (j - 1) ^ 2;24小时(i, j) = r / (i + j - 1);25 H (j,我)= r / (i + j - 1);26日结束27 p = ((n) * p * (n + i)) / (i ^ 2);28日结束

我第一次在机器语言编程算法Burroughs 205数据处理机加州理工学院近60年前。我刚刚从我的青少年。

这是结果n = 6。

格式短T6 = invhilb (6)

T6 = 36 -630 3360 -7560 7560 -2772 -630 14700 -88200 211680 -220500 83160 3360 -88200 564480 -1411200 1512000 -582120 -7560 211680 1512000 -1411200 -3969000 1552320 7560 -220500 1552320 -3969000 4410000 -1746360 -2772 83160 -582120 -1746360 4410000 698544

一个棋盘的信号模式与大整数逆取消希尔伯特矩阵本身的表面光滑。

病人= invhilb (12);间谍(病人> 0)

对数尺度减轻了粗糙度。

冲浪(签署(病人)。*日志(abs(病人)))视图(60 60)

新秀实验

现在使用MATLAB,我要重复这个实验,我做了205年Burroughs当我还是一个新手。我刚刚写了我的第一个程序,用高斯消去法反矩阵。我继续测试它通过反相希尔伯特矩阵和比较结果与精确的逆。(今天的代码使用这个实用工具函数,挑出一个矩阵中最大的元素)。

maxabs = @ (X)马克斯(双(abs (X (:))));

这是n = 10。

n = 10 H = hilb (n);X =发票(H);%计算逆T = invhilb (n);%理论的逆E = X - T;%的错误呃= maxabs (E)

n = 10呃= 5.0259 e + 08年

起初,我可能会说:

哇!错误是10 ^ 8美元。这是一个相当大的误差。我可以相信我的矩阵求逆的代码吗?到底是哪里出了错?

但我知道逆矩阵的元素是巨大的。我们应该看相对错误。

relerr = maxabs (E) / maxabs (T)

relerr = 1.4439 e-04

好的。相对误差是$ 10 ^ {4}$。这似乎仍很多。

我知道希尔伯特矩阵接近奇异。这就是为什么我使用它。约翰·托德是第一个人写条件数据。一个错误估计,包括接近奇点和浮点精度会今天表示

esterr =电导率(H) *每股收益

esterr = 0.0036

这是我理解的。反演过程中的舍入误差放大了矩阵的条件数。,我观察的误差小于估计,这个简单的分析提供了。所以我的反演通过测试的代码。

更深层次的解释

我遇到了吉姆·威尔金森几年后,意识到有更多的故事。实际上我不是反相希尔伯特矩阵。有舍入错误参与计算H之前它是传递到反演程序。

今天,符号数学工具箱可以帮助提供一个更深层次的解释。的“f”国旗在信谊构造函数表示转换为双精度浮点参数完全理性的表现。这是在这种情况下它是如何工作的。为了节省空间,我将打印第一列。

H = hilb (n);F =符号(H,“f”);F (: 1)

ans = 1 1/2 1/4 3602879701896397/18014398509481984 6004799503160661/36028797018963968 2573485501354569/18014398509481984 1/8 2001599834386887/18014398509481984 3602879701896397/36028797018963968 6004799503160661/18014398509481984

的元素H不确切的二进制分数成为大整数的比率。分母是两个的权力;在这种情况下2美元^{54}$和$ 2 ^ {55}$。这些是分子分母除以3美元,5美元等然后四舍五入为最接近的整数。的元素F是接近一个希尔伯特矩阵的具体元素我们可以在二进制浮点数。

让我们反F使用提供的完全理性的算术符号工具箱。(我不能在1960年这么做。)

S =发票(F);

我们现在有三个逆希尔伯特矩阵,X,年代,T。

我们先打印列与对方。

流(% 12 s % 16 s % 15年代\ n”,“X”,“年代”,“T”)流(' % 16.4 f % 16.4 f % 12.0 f \ n”[X (: 1) S (: 1) T (: 1)]”)

X S T 79200 99.9961 -4950 99.9976 100 -4949.6667 - -4949.7926 79192.8929 - 79195.5727 -600535.3362 - -600559.6914 -600600 2522211.3665 -6306300 2522327.5182 2522520 -6305451.1288 - -6305770.4041 9608206.6797 - 9608730.4926 9609600 -8750253.0592 4375800 -8750759.2546 -8751600 4375092.6697 - 4375358.4162 -923624.4113 - -923682.8529 -923780

它看起来像X更接近于年代比年代是T。让我们确认通过计算两个相对误差,之间的区别X和年代和之间的区别年代和T。

格式shorterelerr (1) = maxabs (X - S) / maxabs (T);relerr (2) = maxabs (S - T) / maxabs (T) relerrtotal =总和(relerr)

relerr 9.0252 = 5.4143 e-05 e-05 relerrtotal = 1.4439 e-04

计算中的误差逆来自两个来源,首先生成矩阵,然后计算逆。第一个是大于第二,虽然两相媲美。

得到了MATLAB代码

发表与MATLAB®R2017a