我在 MATLAB 中的实现存在 Frobenius 范数溢出问题答案

【问题标题】：Frobenius norm overflow issue with my implementation in MATLAB我在 MATLAB 中的实现存在 Frobenius 范数溢出问题
【发布时间】：2017-10-13 22:16:33
【问题描述】：

我在 MATLAB 中实现了 Forbenius 范数的以下定义 https://puu.sh/xXsJZ/c42a5e9eac.png

按照描述的方式（即我做了 2 个 for 循环，并通过对每个元素进行平方来增加总和，最后我取总和的平方根）。

我的问题是，有没有办法实现这个规范，这样如果我在矩阵中输入一个相当大的数字，它的平方就不会溢出？在某些情况下，该函数返回“Infinity”，即使真正的 Forbenius 范数可能远低于机器的溢出阈值。（记得在计算结束时取平方根）。

编辑：还有下溢的问题。即使元素 ai,j 不是太小，它的平方也可能下溢。在 MATLAB 中，一个下溢的元素设置为 0。现在只要其他一些元素很大够了，结果还是可以接受的。另一方面，如果所有元素都下溢，我的函数可能会错误地返回 0。

有什么帮助吗？

在我的机器上产生无穷大的矩阵示例是

[1，9999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999999; 1, 1]

但是，当我使用内置的 Frobenius 范数函数时，它可以很好地处理该输入。为什么会这样？

【问题讨论】：

请发布一个重现问题的示例。如果可能，请考虑将这些循环矢量化，您可能会获得很大的速度提升
Matlab 的norm 已经支持 Frobenius 范数。
MATLAB 是如何实现 Frobenius 范数的？我查看了他们的定义并尝试执行 x = sqrt(trace(A* transpose(A)))，但是我仍然遇到溢出问题。当我使用内置 Frobenius 函数 (x = norm(A, 'fro')) 时，我没有遇到溢出问题。

标签： matlab math computer-science

【解决方案1】：

以下代码虽然是用 Java 编写的，但与 Matlab 在计算 Frobenius 范数时所做的很接近。这里的诀窍是hypot 函数，它不仅仅是执行x^2 + y^2，而是计算避免下溢/溢流的斜边。 hypot 在 matlab 中可用，因此不要计算 sqrt(x^2 + y^2)，而是使用 hypot，您应该能够避免下溢/溢流。

public static double normFrob(double[][] matrix) {
    double norm = 0.0;
    int rows = matrix.length;
    int cols = matrix[0].length;
    for (int i = 0; i < rows; ++i) {
        for (int j = 0; j < cols; ++j) {
            norm = hypot(norm, matrix[i][j]);
        }
    }
    return norm;
}

public static double hypot(double a, double b) {
    double r;
    if (Math.abs(a) > Math.abs(b)) {
        r = b / a;
        r = Math.abs(a) * Math.sqrt(1 + r * r);
    } else if (b != 0) {
        r = a / b;
        r = Math.abs(b) * Math.sqrt(1 + r * r);
    } else {
        r = 0.0;
    }
    return r;
}

这段代码 sn-p 取自 Jama，这是为了向 Java 引入一个 Matrix 库，其中一些开发人员是 Matlab 人员。

【讨论】：

Matlab 也有一个内置的hypot 函数。
@horchler 好点。我实际上的目标是伪代码，但碰巧有那个 java sn-p 方便。由于 OP 没有发布任何代码，因此很难说，但我最好的选择是 hypot 将解决他/她的问题。好建议。

【解决方案2】：

溢出保护可以通过多种方式完成。一种可能性是标准化。这是欧几里得范数，但想法是一样的

m_sum = 0;
y=max(abs(x));
for i=1:len(x)
  m_sum = m_sum + (x(i)/y)^2;
end
mynorm = y*sqrt(m_sum);

【讨论】：