希尔伯特矩阵第 12 列后的 Matlab rref() 函数精度误差

Question

我的问题可能很简单，但我想不出对我的问题的合理解释：

当我使用

rref(hilb(8)), rref(hilb(9)), rref(hilb(10)), rref(hilb(11))

它给了我预期的结果，一个单位矩阵。

但是当涉及到

rref(hilb(12))

它没有按预期给出 非奇异 矩阵。我使用了 Wolfram 并且它给出了相同情况下的单位矩阵，所以我确信它应该给出一个单位矩阵。可能存在舍入错误或类似错误，但 1/11 或 1/7 也有一些麻烦的小数

那么为什么 Matlab 在 12 时会这样呢？

Answer 1

这确实像是一个精度错误。这是有道理的，因为 n 阶希尔伯特矩阵的行列式趋于 0，因为 n 趋于无穷大（see here). However, you can use rref with tol parameter：

[R,jb] = rref(A,tol)

并取 tol 非常小以获得更精确的结果。例如，rref(hilb(12),1e-20) 会给你单位矩阵。

编辑-有关tol 参数作用的更多详细信息。

答案底部提供了rref的源代码。 tol用于在某列的某部分搜索绝对值最大的元素，找到主元行。

% Find value and index of largest element in the remainder of column j.
[p,k] = max(abs(A(i:m,j))); k = k+i-1;
   if (p <= tol)
      % The column is negligible, zero it out.
      A(i:m,j) = zeros(m-i+1,1);
      j = j + 1;

如果所有元素的绝对值都小于tol，则该列的相关部分用零填充。这似乎是 rref(hilb(12)) 出现精度错误的地方。通过减少 tol 我们在 rref(hilb(12),1e-20).

中避免了这个问题

源代码：

function [A,jb] = rref(A,tol)
%RREF   Reduced row echelon form.
%   R = RREF(A) produces the reduced row echelon form of A.
%
%   [R,jb] = RREF(A) also returns a vector, jb, so that:
%       r = length(jb) is this algorithm's idea of the rank of A,
%       x(jb) are the bound variables in a linear system, Ax = b,
%       A(:,jb) is a basis for the range of A,
%       R(1:r,jb) is the r-by-r identity matrix.
%
%   [R,jb] = RREF(A,TOL) uses the given tolerance in the rank tests.
%
%   Roundoff errors may cause this algorithm to compute a different
%   value for the rank than RANK, ORTH and NULL.
%
%   Class support for input A:
%      float: double, single
%
%   See also RANK, ORTH, NULL, QR, SVD.

%   Copyright 1984-2005 The MathWorks, Inc. 
%   $Revision: 5.9.4.3 $  $Date: 2006/01/18 21:58:54 $

[m,n] = size(A);

% Does it appear that elements of A are ratios of small integers?
[num, den] = rat(A);
rats = isequal(A,num./den);

% Compute the default tolerance if none was provided.
if (nargin < 2), tol = max(m,n)*eps(class(A))*norm(A,'inf'); end

% Loop over the entire matrix.
i = 1;
j = 1;
jb = [];
while (i <= m) && (j <= n)
   % Find value and index of largest element in the remainder of column j.
   [p,k] = max(abs(A(i:m,j))); k = k+i-1;
   if (p <= tol)
      % The column is negligible, zero it out.
      A(i:m,j) = zeros(m-i+1,1);
      j = j + 1;
   else
      % Remember column index
      jb = [jb j];
      % Swap i-th and k-th rows.
      A([i k],j:n) = A([k i],j:n);
      % Divide the pivot row by the pivot element.
      A(i,j:n) = A(i,j:n)/A(i,j);
      % Subtract multiples of the pivot row from all the other rows.
      for k = [1:i-1 i+1:m]
         A(k,j:n) = A(k,j:n) - A(k,j)*A(i,j:n);
      end
      i = i + 1;
      j = j + 1;
   end
end

% Return "rational" numbers if appropriate.
if rats
    [num,den] = rat(A);
    A=num./den;
end

希尔伯特矩阵第 12 列后的 Matlab rref() 函数精度误差

Matlab rref() function precision error after 12th column of hilbert matrices

precision

matlab