GSoC 2017：实现数值线性代数的迭代求解器

我 GSoC 项目的核心部分是关于在 Julia 中原生实现 Jacobi-Davidson 方法，可在 JacobiDavidson.jl 中找到。此方法计算大型稀疏矩阵 $A$ 和 $B$ 的特征值问题 $Ax = \lambda Bx$ 的一些近似解。由于它在内部使用迭代求解器，因此花费了大量时间来改进 IterativeSolvers.jl。最后，由于迭代求解器通常与预处理器一起使用，因此我还为稀疏矩阵实现了不完全 LU 分解，可在 ILU.jl 中找到。

JacobiDavidson.jl
IterativeSolvers.jl
ILU.jl
示例
1. Jacobi-Davidson
2. ILU 示例
致谢

JacobiDavidson.jl

Jacobi-Davidson 实现已准备好使用，可用于求解非厄米矩阵的（广义）特征值问题。它类似于 Julia 中已有的 eigs 方法：它为您提供几个位于复平面中指定目标附近的特征值。

目前尚未标记任何正式版本，因为仍有一些工作要做：希望广义和普通特征值问题的函数可以在很大程度上合并，因为它们非常相似。此外，还应实现一些针对厄米问题的优化；最后，这些方法尚不支持通用向量和数字。

IterativeSolvers.jl

我们一直在准备 IterativeSolvers.jl 的新版本，该版本提高了 GMRES、CG、切比雪夫迭代、平稳方法和幂方法等求解器的速度和内存使用率。此外，还提供了两种新方法 MINRES 和 BiCGStab(l)，以及针对 Julia 的 SparseMatrixCSC 矩阵类型的平稳方法的高效实现。

此外，该软件包已升级到 Julia 0.6，并且文档结构已重新组织并得到改进。

ILU.jl

线性系统 $Ax = b$ 的迭代方法（例如 BiCGStab(l)）可能不会在任何给定矩阵 $A$ 上快速收敛。通常，如果矩阵 $A$ 只是单位矩阵的扰动，则收敛效果最佳。如果不是这种情况，预处理器可能会有所帮助：与其求解 $Ax = b$ ，您可以尝试求解 $(PA)x = Pb$ ，其中 $P$ 是一个预处理器，使得 $PA$ 更接近单位矩阵。

完美的预处理器将计算 $A$ 的完整 LU 分解，但这需要过多的计算工作，并且需要大量的内存。一个众所周知的技巧是仅近似计算 LU 分解，在过程中丢弃较小的项。这称为不完全 LU 或 ILU。

由于 Julia 中尚未提供 SparseMatrixCSC 类型的 ILU，因此我基于 Na Li、Yousef Saad 和 Edmond Chow 的文章“一般稀疏矩阵的 ILU 的 Crout 版本”对其进行了实现。

软件包 ILU.jl 完全可以使用，并且经过了充分测试。

示例

下面您可以找到一些关于如何使用我一直在开发的软件包的示例。

Jacobi-Davidson

让我们看一个广义特征值问题 $Ax = \lambda Bx$ 的玩具示例，其中 $A$ 和 $B$ 是大小为 $n \times n$ 的对角矩阵，其中 $A_{kk} = \sqrt{k}$

我们使用LinearMaps.jl以无矩阵的方式实现矩阵 $A$ 和 $B$ 的作用。

using LinearMaps

function myA!(y, x)
  for i = 1 : length(x)
    @inbounds y[i] = sqrt(i) * x[i]
  end
end

function myB!(y, x)
  for i = 1 : length(x)
    @inbounds y[i] = x[i] / sqrt(i)
  end
end

n = 100_000
A = LinearMap{Complex128}(myA!, n; ismutating = true)
B = LinearMap{Complex128}(myB!, n; ismutating = true)

矩阵的阶数为100_000。事实证明，如果我们针对频谱内部的特征值，则Jacobi-Davidson内部使用的迭代求解器可能难以解决非常不定的系统。

在这种情况下，我们应该为 $(A - \tau B)$ 使用预处理器，其中 $\tau$ 是目标。我们将只使用精确的逆，它是一个对角矩阵 $P$ ，其元素为 $P_{kk} = \sqrt{k} / (k - \tau)$ 。它可以以无矩阵的方式和就地方式实现。

import Base.LinAlg.A_ldiv_B!

struct SuperPreconditioner{numT <: Number}
    target::numT
end

function A_ldiv_B!(p::SuperPreconditioner, x)
    for i = 1 : length(x)
        @inbounds x[i] *= sqrt(i) / (i - p.target)
    end
end

现在我们使用Near目标调用Jacobi-Davidson并传递预处理器。我们使用GMRES作为迭代求解器，但我们也可以使用BiCGStabl(l)。

using JacobiDavidson

τ = 50_000.1 + 0im
target = Near(τ)
P = SuperPreconditioner(τ)

schur, residuals = jdqz(A, B,
    gmres_solver(n, iterations = 10),
    preconditioner = P,
    target = target,
    pairs = 5,
    ɛ = 1e-9,
    min_dimension = 5,
    max_dimension = 10,
    max_iter = 200,
    verbose = true
)

它收敛到特征值49999、50000、50001、50002和50004。

50004.00000000014 + 3.5749921718300463e-12im
49999.999999986496 - 7.348301591250897e-12im
50001.00000000359 - 1.9761169705101647e-11im
49998.99999999998 - 1.0866253642291695e-10im
50002.00000000171 - 2.3559720511618024e-11im

它还没有检测到50003，但这可能会在稍微增加pairs时发生。由于我们的预处理器，Jacobi-Davidson收敛得非常快。

Residual norm

对于任何给定的问题，构建一个如此好的预处理器并不容易，但通常人们倾向于知道在特定类型的問題中什么效果好。如果没有可用的特定预处理器，您可以始终尝试使用通用预处理器，例如ILU。下一节将对此进行说明。

ILU 示例

作为ILU如何使用的示例，我们生成一个非对称带状矩阵，该矩阵具有通常出现在三维问题有限差分方案中的结构。

n = 64
N = n^3
A = spdiagm((fill(-1.0, n - 1), fill(3.0, n), fill(-2.0, n - 1)), (-1, 0, 1))
Id = speye(n)
A = kron(A, Id) + kron(Id, A)
A = kron(A, Id) + kron(Id, A)
x = ones(N)
b = A * x

矩阵 $A$ 的大小为 $64^3 \times 64^3$ 。我们想要使用例如BiCGStab(2)来解决问题 $Ax = b$ ，但事实证明，当问题规模增大时，收敛速度可能会变慢。一个简单的基准测试表明，将问题求解到合理的容差需要大约2.0秒。

> using BenchmarkTools, IterativeSolvers
> my_x = @btime bicgstabl($A, $b, 2, max_mv_products = 2000);
2.051 s
> norm(b - A * my_x) / norm(b)
1.6967043606691152e-9

现在让我们构建ILU分解。

> using ILU
> LU = crout_ilu(A, τ = 0.1)
> nnz(LU) / nnz(A)
2.1180353639352374

使用上述舍入容差 $\tau$ ，我们的ILU分解仅存储大约两倍于原始矩阵的条目，这是合理的。让我们看看当我们再次对求解器进行基准测试时会发生什么，现在使用ILU作为预条件。

> my_x = @btime bicgstabl($A, $b, 2, Pl = $LU, max_mv_products = 2000);
692.187 ms
> norm(b - A * my_x) / norm(b)
2.133397068536056e-9

它以相同的容差将问题求解速度提高了66%。当然，有一个需要注意的地方，即构建预条件本身也需要时间。

> LU = @btime crout_ilu($A, τ = 0.1);
611.019 ms

因此，总的来说，问题求解速度提高了大约36%。但是，如果我们对同一个矩阵有多个右侧，则可以只构建一次预条件并多次使用它。即使矩阵发生轻微变化，您也可以重用ILU分解。后者正是Jacobi-Davidson中发生的情况。

致谢

我真的很感谢我的GSoC导师Andreas Noack在聊天和视频通话中与我进行的多次讨论。此外，我要感谢Julia社区全体成员给予我的热情欢迎，无论是在线还是在2017年JuliaCon上。