发现PC,Matlab,Maple的极限

我曾经非常自豪地拥有一台Intel P4 2.8G超线程技术的计算机，它的内存刚刚扩容到了1GB，在绝大多数情况下，它都静静地等待着我，等待我去利用它的每一点计算能力，至少我时常这样幻想。回首10年前，那时候主流PC的主频是66M，内存是8M。PC本身的计算能力提高了百倍，物理世界如果这样提高的话，十年前的一层楼就是现在上海最高的建筑物。

但是，无论如何它都只是一台PC（我不得不配上D版软件，Matlab的价格就是USD23000!），世界上的限制永远存在，我们从数学上更是可以发现这些限制——也许说计算复杂度上的天然限制更加合适。我回忆起了今年的线性系统理论考试中的一道题，非常简单问题，给了一个4维线性时不变系统的Jordan标准型，写出能控规范型。实际上就是一个简单的数学变换，假设特征根两两相异的话，就是将系统矩阵A去用一个Vandermonde矩阵以及它的逆做线性变换，其中的Vandermonde矩阵的元素由系统矩阵A的特征多项式组成。但是在考试中你很快会发现，四阶的Vandermonde逆是如此的难求，实际上在考试中绝大多数人都无法通过手算获得正确的答案。考试过后，老师说此题有简便解法，但我这个人一向脑子不太会拐弯，而且比较偏爱IBM的四海一家式解决方案，所以让我们看看，在电脑的帮助下，我们可以最多处理到多少维的线性系统。

现在坐在电脑前，解决这个问题：给定Vandermonde矩阵，求逆，来看看有了现代科技的帮助，4维的系统是多么不值一提。

启动Matlab后依次输入：

% 生成特征值依次为1,2,3,4的Vandermonde矩阵
V=vander(1:4)
% 求逆
Vi=inv(V)
% 验证乘积是否为单位阵
Vi*V

可以看到：

ans =

    1.0000         0   -0.0000   -0.0000
         0    1.0000   -0.0000    0.0000
    0.0000         0    1.0000   -0.0000
         0         0         0    1.0000

虽然格式不是很漂亮，0有的还带个正负，但完全没有问题，前后计算几乎不需要时间！

天性乐观的我马上计算了一下14维的vandermondel矩阵，虽然只是增加了10维，但是马上会遇到数值计算中的舍入误差问题，最终的答案为（为了美观我截去了后面6列类似项）：

ans =

    1.0000    0.0000         0   -0.0000    0.0000    0.0000   -0.0000    0.0000  
   -0.0000    1.0000   -0.0000   -0.0000    0.0000    0.0000   -0.0000   -0.0000  
    0.0000   -0.0000    1.0000   -0.0000   -0.0000    0.0000   -0.0000   -0.0000  
   -0.0000   -0.0000   -0.0000    1.0000   -0.0000   -0.0000   -0.0000   -0.0000  
    0.0005   -0.0001   -0.0000   -0.0000    1.0000   -0.0000   -0.0000   -0.0000  
   -0.0090   -0.0002   -0.0000    0.0000    0.0000    1.0000   -0.0000    0.0000  
    0.0664   -0.0005   -0.0000    0.0000   -0.0000   -0.0000    1.0000    0.0000  
   -0.1250   -0.0059   -0.0018   -0.0003   -0.0000   -0.0000   -0.0000    1.0000  
    0.2500   -0.0586   -0.0103   -0.0007   -0.0001   -0.0000   -0.0000    0.0000  
   -5.2500    0.1563    0.0225    0.0025   -0.0000    0.0000    0.0000   -0.0000  
    8.0000    0.5625    0.0996    0.0081    0.0001   -0.0000    0.0000    0.0000  
  -10.5000   -0.4375   -0.0156   -0.0059    0.0000   -0.0000   -0.0000   -0.0000  
   -1.0000   -0.2188   -0.0234   -0.0002   -0.0006   -0.0000    0.0000    0.0000  
    0.7500    0.0469    0.0151    0.0024    0.0002    0.0000    0.0000    0.0000  

可以看到，实际上高阶vandermondel是非常病态的，通过数值算法见长的Matlab无法承受浮点舍入误差累计，单位阵已经是不成样子了，应该是0的地方出现了8和-10！虽然我有一台超快的PC，计算14阶问题不需要时间，但我还是被浮点精度拌住了，十年来PC的浮点精度可没有一点点的提高，安藤的64位芯片也没有办法。可见，即使你拥有一台超快的计算机，没有合适的软件也白搭;即使你拥有了控制界最令人敬仰的Matlab，号称The language of technical computing也得臣服于IEC 60559:1989, Binary Floating-point Arithmetic for Microprocessor Systems 标准之下。
其实我哪里敢贬低Matlab和IEC(国际电工委员会)，两者都是在学术界，甚至只要是受过高等教育的人都非常敬仰的。但是，中学语文老师教我，欲褒先抑，为了隆重介绍下面的这个软件，我自然就只能说说违心的话了。Maple软件远没有Matlab那么有名，我想在数学软件中，知名度只能排第三（至少在国内次序是Matlab,Mathmatica,Maple）。不过最近我才发现在符号计算中，相对有名的Mathmatica可能比不上Maple，至少Maple是我唯一已知可以直接计算多项式矩阵（Matrix Polynomial Algebra）中左右最大公因子的软件。好了，扯远了，下面还是让我们看看使用Maple能不能突破这个浮点天堑：

同样在Maple中依次输入：

# 在Maple中要显示申明要使用线性代数计算了
with(LinearAlgebra):
# 生成特征值依次为1,2,3,4的Vandermonde矩阵
v:=[$1..14]:
M:=VandermondeMatrix(v);
% 求逆（天哪，为什么就没有一个数学命令标准？每个软件的相同操作都有略微不同的名字）
Mi:=MatrixInverse(M);
% 验证乘积是否为单位阵
Mi.M;

Wow，双击返回的写有14*14矩阵答案框，你可以看到完美的标准14阶单位阵！整个计算也是快的惊人（其实只要计算不超过一秒对于人类都是快的惊人），好奇的我马上试验了一下40阶的Vandermonde，同样没问题，总共计算用时3.35秒，40阶的Vandermonde矩阵如此的病态，条件数已经是10的75次方数量级了。

其实在Matlab中有Symbolic Toolbox符号工具箱，照例也可以解决40阶vandermonde矩阵求逆问题，但是我在Matlab 7.0 SP1环境下就是死机，看来Matlab的符号计算内核Maple 8还有不少问题啊。

实际上，我还试过200阶的Vandermonde矩阵，这基本上是你用PC所能达到的计算极限了，总共计算用时3577秒，基本上一个小时，内存需求486388K，小于512兆内存的PC机可能无法达到这个计算极限。无论如何，这已经是一个天文数字了。实际上，200阶的Vandermonde，矩阵中最大元素为：
200^199=80346902212949513777098104617058130126110149689139641765068800000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000000

或者写成浮点的8.03469E+457，实际上这是一个我们很难通过直观想象的一个巨大的数字，丈量夸克到浩瀚的宇宙只需要10的30次方精度即可完全度量，如果一个原子就是一个宇宙的话，这一数字可以构成的这样的180个循环！

乐观的说，尚未读完大学的学生已经能够勉强处理十年前只有顶尖科学家才能处理的计算了，人类的发展真的是日新月异，悲观的看，我如何能够找到比较实用的数学问题来在这个极限下获得社会的认可呢?

现在我只能说：Never stop thinking——Infineon（天哪，连这样的话也被人家抢先说去了，这年头创新真的很难。）