F# 与数学，第二部分：.NET 实用的 BLAS 和 LAPACK

• 下载源代码 - 4.02 MB

引言

线性代数的应用在深度和广度上都令人惊叹。然而，有三个应用（并非没有重叠）反复出现。

• 求解线性方程组
• 线性最小二乘估计
• 特征值和特征向量计算

在第二部分，我们将探讨 .NET 开发者如何使用 BLAS 和 LAPACK 来求解线性方程组，并在此过程中更深入地了解 BLAS 和 LAPACK 的总体组织结构。

背景

在第一部分中，我们了解到 BLAS 包含了线性代数中最基本的核心例程。在第二部分，我们将在此基础上稍作扩展，加入处理复数矩阵的例程。虽然复数在科学和工程领域很常见，但经济学家通常处理的是实数数据（至少他们是这么希望的）。

然而，第二部分的主要重点是了解 LAPACK 的组织结构。我们在第一部分提到，LAPACK 位于一个更高级别的层次，并在需要时调用 BLAS。LAPACK 函数可以分为三类：

• 驱动程序
• 计算例程
• 辅助例程

如您所知，像求解线性方程组这样的问题需要多个步骤。“驱动程序”负责接受必要的参数；然后，驱动程序会将任务分配给执行这些步骤所需的各个计算例程。辅助例程则提供了一些常见的实用任务。理论上，您甚至可以从不使用驱动程序。您的代码只需要自己调用所有必需的计算例程。而且，我得说，还需要正确调用它们。

随着您在 LAPACK 方面经验的增长，您会注意到有些驱动程序被标记为“专家”。这些驱动程序在处理其任务所需的计算处理方面会做出更复杂的决策。

Using the Code

我们将扩展我们在第一部分开始的 .NET 基础库组织。也就是说，对于每个 BLAS/LAPACK 函数，我们在 LapackAPI.fs 中声明本机函数调用，将托管矩阵与本机实现集成在一个与底层 BLAS/LAPACK 函数同名的文件中，然后添加一个用于与 .NET 世界交互的函数。一个是在 fsLAPACK 中供 F# 开发者使用，另一个是在 clsLapack 中供使用任何 CLS 语言的开发者使用。

添加复数矩阵

在 BLAS/LAPACK 的术语中，“C”开头的函数操作的是单精度复数浮点数值类型，“Z”开头的函数操作的是相应的双精度复数浮点数值类型。我们将添加 ZGEMM 和 ZGEMV，它们是“Z”复数“GE”通用“MM”矩阵-矩阵乘法和它的矩阵-向量乘法姐妹函数“MV”。

这几乎是轻而易举的，因为我们基本上只需要在 F# 代码中将“complex”一词替换为“double”一词。但请记住，这是 FSharp.Math 命名空间中定义的复数类型，稍后我们将不得不考虑 C# 和 VB 程序员使用我们库时可用的选项。.NET Framework 仅在 .NET 4.0 中添加了复数支持。早期版本的开发者必须构建“自己动手”的复数解决方案。

我一直未能找到一种方法让 F# 接受元素类型为 System.Numerics.Complex 的 CLS 数组，尽管看起来应该是可行的。因此，我回到了“暴力”方法，将数组数据复制到 F# 复数数组中。这很丑陋、很慢，但至少目前提供了一个可行的解决方案。不过，这提供了一个机会来说明 Array2D 类中一个非常有用的 static 方法 Array2D.iteri。以下代码将元素类型为 System.Numerics.Complex 的输入数组复制到 F# 复数元素的 Array2D 中。

Array2D.iteri (fun i j _ -> (A.[i,j]<- complex (A_In.[i,j]).Real 
(A_In.[i,j]).Imaginary)) A_In    

F# 的 Array 类有一个“iter”方法，它将一个函数应用于数组的每个元素。iteri 方法与之类似，但当它迭代提供的第二个参数（数组）时，索引 i 和 j 会作为函数的参数提供。请注意下划线，表示可能使用第三个参数。

关于测试脚本的说明

在第一部分中，只有两个函数需要测试。现在有很多了。当我们加载，例如，fsLapack.fs 到交互环境进行测试时，我们必须包含 fsLapack 中所有内容所依赖的每一个文件，即使我们只想测试一个函数。而且，我们必须为每个测试脚本这样做。换句话说，如果我们添加了一个新函数并想创建一个新的测试脚本，我们也需要修改所有以前的测试脚本。出于这个原因，所有的 #load 语句和 "open" 语句都被放在了一个名为 includes.fsx 的单独脚本文件中。fsi 批处理文件在运行感兴趣的特定测试脚本之前会运行它。这样，每次我们创建一个新文件时，我们都在 includes 文件中加载并打开它，就完成了。在 Don Symes 的书中，他说你可以从另一个脚本文件调用像 includes.fsx 这样的脚本文件，但我一直未能让它工作。我采用的方法是在启动 fsi 的批处理文件中加载两个文件。

与第一部分一样，每个测试脚本都配有一个启动它的 fsi 批处理文件。

线性方程组

对于求解线性方程组，我们将再次坚持“一般”情况。Powerpack 单元测试代码提供了一个示例，该示例在很大程度上保留在附件代码中，命名为 fsLinearSystemSolve。此函数使用 DGETRF 和 DGETRS 来完成其任务。DGETRS 在接收 DGETRF 确定的 LU 因子后求解线性方程组。这对于演示很有价值，但 LAPACK 有一个驱动例程 DGESV，可以用单个函数调用来求解实系数线性方程组。源代码中的 fsLinearSystemSolve 函数使用了这种技术，并且本质上与 Powerpack 代码中的示例相同。以下代码使用 DGESV 驱动例程：

 let fsDGESV (A: matrix) (B: vector)  = 
        // we are sticking with the simplest case, B as a vector
        // fortunately, this is also a very common case
        let ipiv = Array.zeroCreate A.NumCols
        let pA = pinM A.Transpose 
        // note: fortran will change the transposed matrix, not A
        // we will not, therefore have access to the factorization
        let pB = pinV B
        let pPivots = pinA ipiv
        let info = dgesv.DGESV_C pA.NativeArray pB.NativeArray pPivots.NativeArray 
        freeM pA; freeV pB; freeA pPivots
        (info, B) 

在 DGESV 中，如果分解的对角线值 U(i,i) 正好为零，则 info 的值将大于零。在这种情况下，矩阵是完全奇异的，计算无法继续。您应该注意到，在浮点表示的世界中，“正好为零”这种情况并不经常发生。一个矩阵在数学上可能是奇异的，但舍入误差会阻止任何对角线元素正好为零。对于 info 参数，“接近”不算。

let DGESV(A : FortranMatrix<double>) (B : FortranMatrix<double>) 
(ipiv : NativeArray<int>) = 
let mutable n = A.NumRows
let mutable nhrs = B.NumCols // should always be one in our work
let mutable lda = A.NumCols
let mutable info = 0
let mutable ldb = B.NumRows 
LAPACK.LapackAPI.DGESV(&&n, &&nhrs, A.Ptr, &&lda, ipiv.Ptr, B.Ptr, &&ldb, &&info);
info

let DGESV_C (A : CMatrix<double>) (B : NativeArray<double>) (ipiv : NativeArray<int>) = 
    // right now we are sticking with vectors on the RHS
    Debug.Assert(A.NumCols = B.Length);
    Debug.Assert(ipiv.Length = A.NumRows);
    DGESV A.NativeTranspose (NativeUtilities.nativeArray_as_FortranMatrix_colvec B) ipiv 

运行代码

批处理文件 RunTestsLinearSystems 找到一组简单的三个未知数三个方程的解。运行批处理文件后，您会注意到该函数返回正确的元组（即 info 和结果向量），但输入矩阵和向量已更改。事实上，您应该会观察到输入向量与函数输出相同。几乎总是如此，底层 Fortran 代码在用于输入的相同数组中返回其结果，这通常不是 .NET 的方式。如果您不希望这样，则必须复制输入数组。此处提供的代码除非必要，否则不会复制数组，因为这对于大型数组来说可能既耗时又耗资源，而且复制的决定应留给使用该库的开发者。处理此问题的最佳方法可能是未来重构的重点。

查看 LU 分解

正如我们刚才提到的，在使用 DGESV 时，LU 分解是在后台完成的，我们无法访问结果。如果我们想查看分解，也许只是为了好玩，我们可以。

在这里，我们发现了 Powerpack 单元测试代码中的另一个小缺陷。如果分解中的对角线元素之一正好为零，DGETRF 会返回大于一的 info 值。Powerpack 代码在这种情况下会抛出异常，如果您使用 LU 分解来求解方阵的线性方程组，这可能是合适的。然而，还有其他原因可以获得 LU 分解。对奇异矩阵或矩形矩阵进行 LU 分解是完全有效的。对于矩形矩阵，几乎可以肯定会存在某些对角线元素为零的情况。这不是错误。我已经注释掉了这个错误并返回了 unit，但这应该在将来进行重构。在提供的代码中，我们将仅限于矩形矩阵的列数多于行数的情况。我们这样做是为了一个重要的简化；我们希望确保 LU 中的 L 始终是方阵。

如何解释 DGETRF 的 LU 分解结果并不立即可见。L 和 U 输出以相加的形式返回到输入矩阵 A 中。假设列数多于行数，L 的对角线元素总是为一，并且完全不表示。因此，要获得 L 和 U，我们必须将 A 分为两部分。L 通过创建一个方阵来获得，取 A 对角线以下的元素，然后将对角线元素设为一。其余元素自然为零。U 从对角线以上的所有内容读取。函数 fsExtractL 和 fsExtractU 已在附件代码中提供，用于分割结果矩阵。乘积 LU 得到一个与输入矩阵等效的矩阵，但行可能已交换位置。枢轴数组的作用是允许我们轻松地重建原始矩阵 A。枢轴数组是一维的，但它指定了在枢轴操作期间交换了哪些行。此信息用于构建表示行交换的置换矩阵。

为此，我们在 fsLapack.fs 中创建了一个辅助函数。我们从一个单位矩阵开始，然后实际执行 Fortran 接收到的枢轴数组指定的行操作。这很容易，因为每行只有一个 1。

let fsPermutationFromPivots (ipivot: int[]) = 
let temp = Matrix.toArray2D (Matrix.identity ipivot.Length)
let pivot_aux i j = 
   if (j = ipivot.[i]-1 && i <> ipivot.[j]-1) then 
      temp.[i,j] <- 1.0
      temp.[j,i] <- 1.0
      temp.[i,i] <- 0.0
      temp.[j,j] <- 0.0

Array2D.iteri (fun i j _ -> (pivot_aux i j)) temp 
Matrix.ofArray2D temp      

现在可以通过从 PLU 获取原始矩阵来测试分解。当然，矩阵乘法是从右到左定义的，所以在 F# 中，乘法将是 P*(L*U)。

即将来到您身边的 F# 编译器...

在下一部分中，我们将研究一般矩阵的特征值和特征向量的计算。

历史

• 2011 年 4 月 24 日：初始发布