【发布时间】：2019-08-27 01:21:59
【问题描述】：
自己尝试在Perl中实现矩阵求逆，发现An Efficient and Simple Algorithm for Matrix Inversion（文章只有两页）。
在我尝试在 Perl 中实现它之后，我发现它不起作用。 我花了 很多 时间试图找出问题所在，所以我得出结论
- 算法不正确
- 我误解了算法
- 我的实现不正确
在展示代码之前，这是一个调试会话，其中包含来自Wikipedia: Inverse Matrix 的示例：
DB<229> $m=[[2,5],[1,3]]
DB<230> x invert($m)
pe[0] == 2
(pivot row 0) 2x2:
2.000 2.500
1.000 3.000
(pivot column 0) 2x2:
2.000 2.500
-0.500 3.000
(rest 0) 2x2:
2.000 2.500
-0.500 1.750
(pivot 0) 2x2:
0.500 2.500
-0.500 1.750
pe[1] == 1.75
(pivot row 1) 2x2:
0.500 2.500
-0.286 1.750
(pivot column 1) 2x2:
0.500 -1.429
-0.286 1.750
(rest 1) 2x2:
0.908 -1.429
-0.286 1.750
(pivot 1) 2x2:
0.908 -1.429
-0.286 0.571
0 1
1 3.5
DB<231>
这是我写的代码：
#!/usr/bin/perl -w
use 5.026;
use strict;
# invert matrix
# An Efficient and Simple Algorithm for Matrix Inversion
# Ahmad Farooq, King Khalid University, Saudi Arabia
# Khan Hamid, National University of Computer and Emerging Sciences (NUCES),
# Pakistan
sub invert($)
{
my $m = shift; # matrix is an array of rows
my ($pp, $det);
my ($rp, $pe);
my $n = scalar(@$m);
for ($pp = 0, $det = 1.0; $pp < $n; ++$pp) {
$rp = $m->[$pp]; # pivot row
$pe = $rp->[$pp]; # pivot element
print "pe[$pp] == $pe\n";
last if ($pe == 0); # Epsilon test?
$det *= $pe;
# calculate pivot row
for (my $j = 0; $j < $n; ++$j) {
next if ($j == $pp);
$rp->[$j] /= $pe;
}
pm($m, "pivot row $pp");
# calculate pivot column
for (my $i = 0; $i < $n; ++$i) {
next if ($i == $pp);
$m->[$i]->[$pp] /= -$pe;
}
pm($m, "pivot column $pp");
for (my $j = 0; $j < $n; ++$j) {
next if ($j == $pp);
for (my ($i, $rj) = (0, $m->[$j]); $i < $n; ++$i) {
next if ($i == $pp);
$rj->[$i] += $rp->[$j] * $m->[$i]->[$pp];
}
}
pm($m, "rest $pp");
$rp->[$pp] = 1.0 / $pe;
pm($m, "pivot $pp");
}
return ($pe != 0.0, $det);
}
pm() 函数只是一个用于调试目的的“打印矩阵”：
# print matrix
sub pm($;$)
{
my ($m, $label) = @_;
my $n = scalar(@$m);
print "($label) " if ($label);
print "${n}x${n}:\n";
for (my $i = 0; $i < $n; ++$i) {
for (my $j = 0; $j < $n; ++$j) {
if (defined(my $v = $m->[$i]->[$j])) {
printf('%8.3f', $v);
} else {
print ' ???????';
}
}
print "\n";
}
}
有什么见解吗？
复制提示（添加于 2019-08-28）
我认为这很明显，但以防万一： 如果你想重现调试会话中显示的输出，也许只需在代码末尾添加这两行：
my $m=[[2,5],[1,3]]; # matrix to invert
print join(', ', invert($m)), "\n"; # invert $m, printing result
注意（添加于 2019-09-02）：
对于 Wikipedia 文章 ($m = [[1, 2, 0], [2, 4, 1], [2, 1, 0]]) 中给出的 3x3 矩阵，该算法失败，因此真正的实现应该转向改进的算法（可以选择对角线之外的枢轴元素）。
【问题讨论】：
-
use strict在use 5.026之后是多余的，但是您缺少use warnings。 -
也许您可以将示例输入矩阵和反转矩阵与脚本的输出一起发布，以便我们查看是否存在错误模式。
-
@lordadmira：你能说出第一个代码块中给出的例子有什么问题吗？
-
抱歉，我并没有真正关注调试器的输出。不幸的是，我还没有时间研究算法。
-
我认为问题出在参考论文的第 7 步。请注意，
a[i,p]上有一个素数，但a[p,j]上没有。所以我猜你必须保存a[p,j]的元素？