标签: numerical-stability

我正在尝试为Python编写Wilson的谱密度分解算法[1]的实现.该算法迭代地将[QxQ]矩阵函数分解为其平方根(它是用于谱密度矩阵的Newton-Raphson平方根查找器的扩展).

问题是我的实现仅收敛于45x45及更小的矩阵.因此,在20次迭代之后,矩阵之间的总平方差约为2.45e-13.但是,如果我输入大小为46x46的输入,则直到第100次迭代才会收敛.对于47x47或更大,矩阵永远不会收敛; 对于大约100次迭代,误差在100和1000之间波动,然后开始非常快速地增长.

您将如何尝试调试此类内容？似乎没有任何特定的点,它变得疯狂,并且矩阵太大,我实际上不能手动进行计算.有没有人有寻找奇怪的数字错误的提示/教程/启发式？

我以前从未处理过这样的事情,但我希望你们有些人......

谢谢, - 丹

[1] GT Wilson."基质谱密度的因式分解".SIAM J. Appl.数学(Vol 23,No.4,1972年12月)

我一直在试图计算只有三个不同的项目,组成一个36分36矩阵的Jordan标准型1,1/2和0.矩阵是概率转移矩阵,因此,给定这些条目,矩阵显然是稀疏的.

我一直遇到的问题如下:每当我尝试计算时

[V, J] = jordan(A),

要么

[V, J] = jordan(sym(A)),

我收到以下错误消息:

在MuPAD命令中使用mupadmex 错误时出错:相似度矩阵太大.

sym/mupadmexnout出错(第1546行)
        out = mupadmex(fcn,args {:});

sym/jordan中的错误(第32行)
        [Vsym,Jsym] = mupadmexnout('symobj :: jordan',A,'All');

我在MATLAB帮助中读到,Jordan形式的计算对扰动非常敏感.但是,我认为我的计算不是问题,因为矩阵的所有条目都是整数或整数比.

我的问题如下:

1. 如何解释收到的错误输出？
2. 我收到的错误是否可以解决？
3. 如果错误不可寻址,是否有替代方法(Matlab中的函数)我可以尝试计算Jordan形式？

我需要在Matlab中计算向量的归一化指数.

简单写作

res = exp(V)/sum(exp(V))

V元素中的溢出大于log(realmax)= 709.7827.(我不确定下溢情况.)

我该如何实现它以避免数值不稳定？

更新:我收到了关于如何避免溢出的出色回复.但是,我仍然很高兴听到您对代码中可能出现下溢的想法.

我猜有一些我无法找到的标准技巧:无论如何,我想以数字稳定的方式计算一个非常接近1的数字的大功率(认为1-p,其中p <1e-17) .在我的系统上,1-p被截断为1.

使用对数的泰勒展开式,我得到以下边界

我能做些什么更聪明的事吗？

添加是否x + x可以通过IEEE 754(IEC 559)浮点标准中的乘法2 * x进行互换,或者更一般地说,是否有任何保证并始终给出完全相同的结果？case_addcase_mul

#include <limits>

template <typename T>
T case_add(T x, size_t n)
{
    static_assert(std::numeric_limits<T>::is_iec559, "invalid type");

    T result(x);

    for (size_t i = 1; i < n; ++i)
    {
        result += x;
    }

    return result;
}

template <typename T>
T case_mul(T x, size_t n)
{
    static_assert(std::numeric_limits<T>::is_iec559, "invalid type");

    return x * static_cast<T>(n);
}

如果我有一个随机数Z，它被定义为另外两个随机数 和 的总和，X那么Y的概率分布是和Z的概率分布的卷积。卷积基本上是分布函数乘积的积分。卷积中的积分通常没有解析解，因此必须使用基本求积算法来计算。在伪代码中：XY

prob_z(z) = integrate(lambda t: prob_x(t) * prob_y(z-t), -inf, inf)\n

举一个具体的例子，可以使用以下Python/Scipy代码计算Z正态分布变量X和对数正态分布变量的总和：Y

from scipy.integrate import quad\nfrom scipy.stats import norm, lognorm\nfrom scipy import log\n\nprob_x = lambda x: norm.pdf(x, 0, 1)  # N(mu=0, sigma=1)\nprob_y = lambda y: lognorm.pdf(y, 0.1, scale=10)  # LogN(mu=log(10), sigma=0.1)\ndef prob_z(z):\n    return quad(lambda t: prob_x(t)*prob_y(z-t), -inf, inf)\n

现在我想计算对数概率。天真的解决方案是简单地执行以下操作：

def log_prob_z(z):\n    return log(prob_z(z))\n

然而，这在数值上是不稳定的。在大约 39 个标准差之后，概率分布在数值上为 …

我正在使用 CVXPY（版本 1.0）来解决二次程序（QP），我经常遇到这个异常：

求解器错误：求解器“xxx”失败。尝试另一个求解器。

这使我的程序非常脆弱。我尝试了不同的求解器，包括 CVXOPT、OSQP、ECOS、ECOS_BB、SCS。他们都或多或少有相同的问题。我注意到，当我使求解器的停止标准更严格（例如，降低绝对误差容限）时，我会变SolverError得更频繁，而当我使其不那么严格时，SolverError问题会减弱甚至消失。我还发现 CVXPY 抛出的方式SolverError是随机的：如果我多次运行同一个程序，有一些运行SolverError会获得最佳结果而另则会获得最佳结果。

虽然我可以通过尝试更多次并降低停止标准来避免 SolverError，但我真的很想了解异常背后的真正具体原因

求解器错误：求解器“xxx”失败。尝试另一个求解器。

这个错误并没有真正提供信息，我不知道如何提高解决问题的稳健性。其原因是否特定于求解器？是否为一组明确定义的情况抛出此异常？或者它只是一种说“由于未知原因出现问题”的方式？这些可能是什么原因？

我知道如何通过添加到元素 -max _i x_i 来使 softmax 稳定。这避免了上溢和下溢。现在，记录此日志可能会导致下溢。log softmax(x) 可以计算为零，导致 -infinity。

我不知道如何解决它。我知道这是一个常见问题。我读了几个答案，我不明白。但是我仍然对如何解决这个问题感到困惑。

PS：如果你提供一个简单的例子，那就太棒了。

我正在尝试使用 Python 研究以下延迟微分方程的行为：

y\'\'(t) = -y(t)/\xcf\x84^2 - 2y\'(t)/\xcf\x84 - Nd*f(y(t-T))/\xcf\x84^2,\n

其中f是截止函数，当其参数的绝对值在 1 到 10 之间时，它本质上等于恒等式，否则等于 0（见图 1），并且 、 和Nd是\xcf\x84常数T。

为此，我使用 JiTCDDE 包。这为上式提供了合理的解。尽管如此，当我尝试在方程右侧添加噪声时，我得到了一个在几次振荡后稳定为非零常数的解。这不是方程的数学解（唯一可能的常数解等于零）。我不明白为什么会出现这个问题以及是否可以解决它。

我在下面重现我的代码。这里，为了简单起见，我用高频余弦代替噪声，将其引入方程组作为虚拟变量的初始条件（余弦可以直接引入系统，但对于一般噪音这似乎不可能）。为了进一步简化问题，我还删除了涉及该函数的术语f，因为没有它也会出现问题。图 2 显示了代码给出的函数图。

from jitcdde import jitcdde, y, t\nimport numpy as np\nfrom matplotlib import pyplot as plt\nimport math\nfrom chspy import CubicHermiteSpline\n\n\n# Definition of function f:\ndef functionf(x):\n    return x/4*(1+symengine.erf(x**2-Bmin**2))*(1-symengine.erf(x**2-Bmax**2))\n\n#parameters:\n\xcf\x84 = 42.9\nT = 35.33\nNd = 8.32\n\n# Definition of the initial conditions:\ndt = .01 …

我有一个Tensorflow多类分类器，其产生nan或inf同时使用计算概率tf.nn.softmax。请参阅以下代码段（logits具有形状batch_size x 6，因为我有 6 个类并且输出是单热编码的）。batch_size是 1024。

logits = tf.debugging.check_numerics(logits, message='bad logits', name=None)
probabilities = tf.nn.softmax(logits=logits, name='Softmax')
probabilities = tf.debugging.check_numerics(probabilities, message='bad probabilities', name=None)

分类器在找到nan或inf中的最后一条语句时失败probabilities。logits是干净的，否则第一个语句将失败。

从我读到的内容来看tf.nn.softmax，它可以处理 logits 中非常大和非常小的值。我已经在交互模式下验证了这一点。

>>> with tf.Session() as s:
...   a = tf.constant([[1000, 10], [-100, -200], [3, 4.0]])
...   sm = tf.nn.softmax(logits=a, name='Softmax')
...   print(a.eval())
...   print(sm.eval())
...
[[1000.   10.]
 [-100. -200.]
 [   3.    4.]]
[[1. …