使浮点数低于某个公差 0.0

Question

经过一些计算，我有一个浮点数矩阵，其中有一些非常小的 negative 条目（例如 -2.6e-9），它们近似于 0。

我知道这个矩阵应该是正定的，即所有正特征值，但是，我怀疑这些大约为 0 但负条目使矩阵不是正定的，因此我不能做我想做的进一步计算（Cholesky 分解，对于那些想知道的人）。

我的问题是这个：无论如何我可以强制这些非常小的条目负条目，它们近似于 0 为 0.0？

我想我在这里误解了浮点数！

Answer 1

不一定是最快的，但似乎可行：

julia> tol = 1e-5
1.0e-5

julia> X = [0.1 -1e-4 -1e-20; -2.6e-9 0.1 1.1e5]
2×3 Matrix{Float64}:
  0.1     -0.0001      -1.0e-20
 -2.6e-9   0.1     110000.0

julia> Y = ifelse.((X .< -tol ) .|  (X .>  0.0 ), X, 0.0)
2×3 Matrix{Float64}:
 0.1  -0.0001       0.0
 0.0   0.1     110000.0

这样做的好处是保持真正的负数为负数，因此您可以发现最终的错误。

Answer 2

A[A .≈ 0.0] .= 0.0

≈，可以用latex快捷键输入\approx<tab>，shorthand是isapprox函数，可以直接用它来设置具体的公差希望改用。

编辑：如前所述，默认值不适用于针对 0.0 的测试。

A[isapprox.(A, 0.0, atol = 1e-8)] .= 0.0

Answer 3

当然，使用 isapprox 可以实现具有任意容差的一种相对简单的方法。例如：

julia> A = fill(-2.6e-9, 5,5)
5×5 Matrix{Float64}:
 -2.6e-9  -2.6e-9  -2.6e-9  -2.6e-9  -2.6e-9
 -2.6e-9  -2.6e-9  -2.6e-9  -2.6e-9  -2.6e-9
 -2.6e-9  -2.6e-9  -2.6e-9  -2.6e-9  -2.6e-9
 -2.6e-9  -2.6e-9  -2.6e-9  -2.6e-9  -2.6e-9
 -2.6e-9  -2.6e-9  -2.6e-9  -2.6e-9  -2.6e-9

julia> map!(x -> isapprox(x, 0, atol=1e-8) ? 0 : x, A, A)
5×5 Matrix{Float64}:
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0
 0.0  0.0  0.0  0.0  0.0

请注意，使用 map! 执行此操作不会分配，因此应该非常有效：

julia> @btime map!(x -> isapprox(x, 0, atol=1e-8) ? 0 : x, $A, $A)
  30.984 ns (0 allocations: 0 bytes)

如果您不指定公差（上面 atol 选项的作用），它将使用稍微复杂的公式为一般浮点运算选择一个合理的默认值：

  isapprox(x, y; atol::Real=0, rtol::Real=atol>0 ? 0 : √eps, nans::Bool=false[, norm::Function])

  Inexact equality comparison. Two numbers compare equal if their relative distance or
  their absolute distance is within tolerance bounds: isapprox returns true if norm(x-y)
  <= max(atol, rtol*max(norm(x), norm(y))). The default atol is zero and the default
  rtol depends on the types of x and y. The keyword argument nans determines whether or
  not NaN values are considered equal (defaults to false).

  For real or complex floating-point values, if an atol > 0 is not specified, rtol
  defaults to the square root of eps of the type of x or y, whichever is bigger (least
  precise). This corresponds to requiring equality of about half of the significand
  digits. Otherwise, e.g. for integer arguments or if an atol > 0 is supplied, rtol
  defaults to zero.

关于您之前的计算是否有问题：不一定，因为舍入误差是使用固定精度浮点数不可避免的部分。然而 -2.6e-9 与

相比确实听起来有点大

julia> eps(0.0)
5.0e-324

所以这取决于你到底在做什么。例如，如果您要从彼此中减去两个大但相似的浮点数，那么对于 Float64s

这种错误级别可能并不令人惊讶

Answer 4

你的问题有几个很好的答案，所以让我评论一下如果你有舍入误差如何使你的矩阵正定。最简单的方法是在原始矩阵中添加一个小的对角矩阵。让我举个例子。采用以下矩阵：

julia> A = fill(-2.6e-12, 5, 5)
5×5 Matrix{Float64}:
 -2.6e-12  -2.6e-12  -2.6e-12  -2.6e-12  -2.6e-12
 -2.6e-12  -2.6e-12  -2.6e-12  -2.6e-12  -2.6e-12
 -2.6e-12  -2.6e-12  -2.6e-12  -2.6e-12  -2.6e-12
 -2.6e-12  -2.6e-12  -2.6e-12  -2.6e-12  -2.6e-12
 -2.6e-12  -2.6e-12  -2.6e-12  -2.6e-12  -2.6e-12

您不会通过使所有小的负值等于 0 来使其成为正定的。但是，如果你在它的对角线上添加一些小值，它就会变成正定的：

julia> A + sqrt(eps())*I
5×5 Matrix{Float64}:
  1.48986e-8  -2.6e-12     -2.6e-12     -2.6e-12     -2.6e-12
 -2.6e-12      1.48986e-8  -2.6e-12     -2.6e-12     -2.6e-12
 -2.6e-12     -2.6e-12      1.48986e-8  -2.6e-12     -2.6e-12
 -2.6e-12     -2.6e-12     -2.6e-12      1.48986e-8  -2.6e-12
 -2.6e-12     -2.6e-12     -2.6e-12     -2.6e-12      1.48986e-8

julia> isposdef(A + sqrt(eps())*I)
true

在实践中，您可能需要做两件事：将小的负值四舍五入为 0（如果您知道在您的问题中不可能输入负数）并在您的对角线上添加一些小值矩阵.

Answer 5

有一个包可以在矩阵上强制正定性，称为 Positive Factorizations.jl。

鉴于您有一个 H 矩阵，您可以通过以下方式进行正 Cholesky 分解：

using PositiveFactorizations
A = fill(-2.6e-12, 5, 5)
F = cholesky(Positive, A)

使浮点数低于某个公差 0.0

Making floating point numbers below some tolerance 0.0

floating-point

precision

numbers

julia