调试线性代数相关问题#
线性代数相关的错误报告是诊断和解决中最具挑战性的问题之一。这不仅是因为线性代数在数学/算法上可能具有挑战性（这对 SciPy 的许多部分来说都是如此），而且还因为 BLAS/LAPACK 库是一个复杂的构建时以及运行时依赖项 - 而且我们支持大量的 BLAS/LAPACK 库。
本文档旨在提供有关如何调试线性代数问题的指导。
如果存在真正的错误，它可能存在于以下三个位置之一
在正在使用的 BLAS 库中，
在 SciPy 的 BLAS 或 LAPACK 绑定中（由
numpy.f2py和/或 Cython 生成），在 SciPy 的算法代码中。
第一个关键步骤是确定错误是在 SciPy 中还是在 BLAS 库中。为此，最有效的方法是设置您的环境，使其能够在不同的 BLAS 库之间实现运行时切换（这是我们不提供开箱即用支持的，并且使用 PyPI 的 SciPy wheels 是不可能的）。
上游 BLAS 库的作者强烈希望获得干净的复现器（就像我们一样），这意味着：不涉及 Python。因此，本指南还将介绍如何在 C 或 Fortran 中创建复现器。
查找正在使用的 BLAS 库#
SciPy 有一个函数，show_config，用于自省已安装软件包的构建配置。这允许查询 BLAS 和 LAPACK 的详细信息。例如
>>> blas_dep = scipy.show_config(mode='dicts')['Build Dependencies']['blas']
>>> for key in blas_dep:
... print(f"{key}: {blas_dep[key]}")
...
name: openblas
found: True
version: 0.3.23
detection method: pkgconfig
include directory: /home/user/miniforge/envs/scipy-dev/include
lib directory: /home/user/miniforge/envs/scipy-dev/lib
openblas configuration: USE_64BITINT=0 DYNAMIC_ARCH=1 DYNAMIC_OLDER= NO_CBLAS= NO_LAPACK=0 NO_LAPACKE= NO_AFFINITY=1 USE_OPENMP=0 PRESCOTT MAX_THREADS=128
pc file directory: /home/user/miniforge/envs/scipy-dev/lib/pkgconfig
此方法对于 SciPy wheels 和本地开发构建是正确的。它可能对其他安装是正确的，但是请记住，像 conda-forge 和 Debian 这样的发行版可能会针对存根库（通常是 blas.so/lapack.so）进行构建，然后为用户安装另一个库 - 在这种情况下，即使例如在环境中安装了 OpenBLAS 或 MKL，也会报告普通的 blas 和 lapack。对于此类安装，threadpoolctl 通常能够报告正在使用的实际 BLAS 库（除了它不报告普通的 Netlib BLAS，请参阅 threadpoolctl#159）
$ python -m threadpoolctl -i scipy.linalg
[
{
"user_api": "blas",
"internal_api": "openblas",
"prefix": "libopenblas",
"filepath": "/home/user/miniforge/envs/dev/lib/libopenblasp-r0.3.21.so",
"version": "0.3.21",
"threading_layer": "pthreads",
"architecture": "SkylakeX",
"num_threads": 24
}
]
在本地开发环境中可能有帮助的其他自省方法包括
检查共享库的依赖项
$ ldd build/scipy/linalg/_fblas.cpython-*.so
...
libopenblas.so.0 => /home/user/miniforge/envs/scipy-dev/lib/libopenblas.so.0 (0x0000780d6d000000)
% otool -L build/scipy/linalg/_fblas.cpython-310-darwin.so
build/scipy/linalg/_fblas.*.so:
/usr/lib/libSystem.B.dylib (compatibility version 1.0.0, current version 1336.61.1)
@rpath/libopenblas.0.dylib (compatibility version 0.0.0, current version 0.0.0)
检查链接的库是否包含给定的符号。例如，conda-forge 安装一个
libblas.so，它可以是任何受支持的库
$ nm -gD ~/miniforge/envs/scipy-dev/lib/libblas.so | rg openblas_set_num_threads
0000000000362990 T openblas_set_num_threads
% nm ~/miniforge/envs/scipy-dev/lib/libblas.3.dylib | rg openblas_set_num_threads
000000000015b6b0 T _openblas_set_num_threads
设置环境以切换 BLAS 库#
我们将介绍几种在不同 BLAS 库之间切换的方法，因为最简单的方法可能取决于您的操作系统/发行版以及您想要 SciPy 的发布版本还是开发版本。
Conda-forge#
也许最简单的方法是使用发行版提供的运行时切换功能。例如，要创建一个安装了最新 SciPy 版本的 conda 环境，然后在 OpenBLAS、Netlib BLAS/LAPACK 和 MKL 之间切换，就像
$ mamba create -n blas-switch scipy threadpoolctl
$ mamba activate blas-switch
$ python -m threadpoolctl -i scipy.linalg
...
"user_api": "blas",
"internal_api": "openblas",
$ mamba install "libblas=*=*netlib"
...
libblas 3.9.0-21_linux64_openblas --> 3.9.0-5_h92ddd45_netlib
libcblas 3.9.0-21_linux64_openblas --> 3.9.0-5_h92ddd45_netlib
liblapack 3.9.0-21_linux64_openblas --> 3.9.0-5_h92ddd45_netlib
$ mamba install "libblas=*=*mkl"
...
libblas 3.9.0-5_h92ddd45_netlib --> 3.9.0-21_linux64_mkl
libcblas 3.9.0-5_h92ddd45_netlib --> 3.9.0-21_linux64_mkl
liblapack 3.9.0-5_h92ddd45_netlib --> 3.9.0-21_linux64_mkl
$ python -m threadpoolctl -i scipy.linalg
...
"user_api": "blas",
"internal_api": "mkl",
当通过 dev.py 构建时，也可以对开发构建执行此操作，方式与 SciPy 的 conda-forge 构建配方 完全相同（为简洁起见，省略了输出，它们与上面的输出类似）
$ mamba env create -f environment.yml
$ mamba activate scipy-dev
$ mamba install "libblas=*=*netlib" # necessary, we need to build against blas/lapack
$ python dev.py build -C-Dblas=blas -C-Dlapack=lapack -C-Duse-g77-abi=true
$ python dev.py test -s linalg # run tests to verify
$ mamba install "libblas=*=*mkl"
$ python dev.py test -s linalg
$ mamba install "libblas=*=*openblas"
Linux 发行版包管理器#
许多 Linux 发行版使用 update-alternatives 机制，允许通过系统包管理器在不同的 BLAS 库之间切换。请注意，这是一种通用的机制，用于管理“同一库或工具的多个实现”情况，而不是 BLAS/LAPACK 特有的机制。它类似于上面的 conda-forge 方法，因为它适用于发行版提供的 scipy 包以及针对引用 libblas/liblapack 接口的开发构建。
接口如下所示
$ update-alternatives --config libblas.so.3
$ update-alternatives --config liblapack.so.3
这将在终端中打开一个菜单，其中包含所有可供选择的可用库。由于接口和可用选项可能因发行版而异，因此我们在此链接到 Debian 关于 BLAS/LAPACK 切换的文档，避免更详细地记录此功能在其他发行版上的工作方式。
请注意，Fedora 是一个例外；它是唯一一个发布 FlexiBLAS（有关详细信息，请参阅下一节）的发行版，并允许并行安装多个 BLAS 库，因此无需调用系统包管理器即可实现真正的运行时切换。有关更多详细信息，请参阅 Fedora 关于系统级和用户级选择的文档。
FlexiBLAS#
FlexiBLAS 为它可以检测到的所有已安装的 BLAS 库提供运行时切换支持（以及其他功能）。在撰写本文时（2024 年 3 月），存在一些限制，主要是：不支持 Windows，不支持 macOS Accelerate（更新的版本，带有 NEWLAPACK 符号）。如果这些限制对您无关紧要，则 FlexiBLAS 可以成为一个非常有用的高效调试工具，包括对于您必须从源代码构建的 OpenBLAS 和其他 BLAS 库的版本。
一旦您设置好一切，开发体验将是
$ python dev.py build -C-Dblas=flexiblas -C-Dlapack=flexiblas
$ FLEXIBLAS=NETLIB python dev.py test -s linalg
$ FLEXIBLAS=OpenBLAS python dev.py test -s linalg
# Or export the environment variable to make the selection stick:
$ export FLEXIBLAS=OpenBLAS
您还可以提供到构建的 BLAS 库的路径（例如，FLEXIBLAS="libbhlas_atlas.so"）- 有关更多详细信息，请参阅 其 README 中的用法文档。
除非您使用的是 Fedora，否则您可能必须从源代码构建 FlexiBLAS，这需要做一些工作。好消息是，无论您使用的是 Linux 还是 macOS，以及通过 virtualenvs、conda 环境还是以其他方式使用 Python，这都应该有效。我们将介绍如何从源代码构建 OpenBLAS 和 FlexiBLAS，以便调试最新 OpenBLAS 版本中的某些内容是否与 Netlib BLAS/LAPACK（或 MKL）不同。
以下内容应在您可以构建 SciPy 本身的任何环境中工作；我们唯一需要的额外工具是 CMake（例如，使用 pip install cmake 安装）。
克隆每个存储库
$ cd .. # starting from the root of the local `scipy` repo
$ mkdir flexiblas-setup && cd flexiblas-setup
$ git clone https://github.com/OpenMathLib/OpenBLAS.git openblas
$ git clone https://github.com/mpimd-csc/flexiblas.git
$ mkdir built-libs # our local prefix to install everything to
构建 OpenBLAS
$ cd openblas
$ mkdir build && cd build
$ cmake .. -DBUILD_SHARED_LIBS=ON -DCMAKE_INSTALL_PREFIX=$PWD/../../built-libs
$ cmake --build . -j
$ cmake --install . --prefix $PWD/../../built-libs
$ cd ../..
构建 FlexiBLAS
$ cd flexiblas
$ mkdir build && cd build
$ # Note: this will also pick up the libraries in your system/env libdir
$ cmake .. -DEXTRA="OpenBLAS" -DLAPACK_API_VERSION=3.9.0 \
-DOpenBLAS_LIBRARY=$PWD/../../built-libs/lib/libopenblas.so \
-DCMAKE_INSTALL_PREFIX=$PWD/../../built-libs
$ cmake --build . -j
$ cmake --install . --prefix $PWD/../../built-libs
$ cd ../..
我们现在准备好针对 FlexiBLAS 构建 SciPy
$ export PKG_CONFIG_PATH=$PWD/flexiblas-setup/built-libs/lib/pkgconfig/
$ cd scipy
$ python dev.py build -C-Dblas=flexiblas -C-Dlapack=flexiblas
...
Run-time dependency flexiblas found: YES 3.4.2
现在我们可以运行测试。请注意，构建 NETLIB 选项时无需指定；它是 FlexiBLAS 中的默认值，并且源代码包含在其存储库中
$ FLEXIBLAS=OpenBLAS python dev.py test -s linalg
$ FLEXIBLAS=NETLIB python dev.py test -s linalg
$ python dev.py test -s linalg # uses the default (NETLIB)
此后端切换也可以在 Python 解释器中使用 threadpoolctl 完成（有关详细信息，请参阅 其 README）。
可以使用以下命令检查系统上可用的其他库
$ ./flexiblas-setup/built-libs/bin/flexiblas list
注意
使用参考 BLAS/LAPACK 或 BLIS 的本地构建更加困难，因为 FlexiBLAS 需要一个包含所有必要符号的单个共享库。将单独的 libblas 和 liblapack 用作“系统库”可能是可行的，但这已被证明更脆弱且难以构建（因此这是 YMMV）。如果您确实想尝试
构建参考 BLAS 和 LAPACK
$ git clone Reference-LAPACK/lapack.git $ cd lapack $ mkdir build && cd build $ cmake -DCBLAS=ON -DBUILD_SHARED_LIBS=OFF .. $ cmake –build . -j $ cmake –install . –prefix $PWD/../../built-libs
然后将以下两行添加到 FlexiBLAS 的 cmake .. 配置命令中
-DSYS_BLAS_LIBRARY=$PWD/../../built-libs/lib/libblas.a \
-DSYS_LAPACK_LIBRARY=$PWD/../../built-libs/lib/liblapack.a \
在 C 或 Fortran 中创建复现器#
我们的经验告诉我们，绝大多数错误都存在于 SciPy 内部，而不是在 OpenBLAS 或其他 BLAS 库中。但是，如果使用不同的 BLAS 库进行的测试告诉我们问题是特定于单个 BLAS 库的（甚至可能是该库的单个版本存在回归），那么下一步是在 C 或 Fortran 中生成复现器；这样做对于向上游报告错误是必要的，并且使 BLAS 库开发人员更容易解决问题。
要从使用 NumPy 数组作为输入的 Python 复制器（使用 scipy 函数）转换为 C/Fortran 复制器，需要找到 SciPy 中执行的代码路径，并确定调用了哪个确切的 BLAS 或 LAPACK 函数，以及使用了什么输入（注意：答案可能包含在 .pyf.src f2py 签名文件中；查看生成目录中的 _flapackmodule.c 也可能很有用）。然后，可以通过定义一些整数/浮点变量和数组（通常，使用硬编码数字的小数组就足够了）在 C/Fortran 中重现此过程。
可以在例如 Netlib LAPACK 文档或 Reference-LAPACK/lapack 存储库中查找 BLAS 和 LAPACK 函数的参数列表。
下面显示了一个参考 LAPACK 中问题的复制器，该问题在 scipy#11577 中被报告为一个 SciPy 问题。我们将文件命名为 ggev_repro_gh_11577.c|f90
#include <stdio.h>
#include "lapacke.h"
#define n 4
int main()
{
int lda=n, ldb=n, ldvr=n, ldvl=n, info;
char jobvl='V', jobvr='V';
double alphar[n], alphai[n], beta[n];
double vl[n*n], vr[n*n];
// int lwork = 156;
// double work[156]; /* cheat: 156 is the optimal lwork from the actual lapack call*/
double a[n*n] = {12.0, 28.0, 76.0, 220.0,
16.0, 32.0, 80.0, 224.0,
24.0, 40.0, 88.0, 232.0,
40.0, 56.0, 104.0, 248.0};
double b[n*n] = {2.0, 4.0, 10.0, 28.0,
3.0, 5.0, 11.0, 29.0,
5.0, 7.0, 13.0, 31.0,
9.0, 11.0, 17.0, 35.0};
info = LAPACKE_dggev(LAPACK_ROW_MAJOR, jobvl, jobvr, n, a, lda, b,
ldb, alphar, alphai, beta, vl, ldvl, vr, ldvr); //, work, lwork, info);
printf("info = %d\n", info);
printf("Re(eigv) = ");
for(int i=0; i < n; i++){
printf("%f , ", alphar[i] / beta[i] );
}
printf("\nIm(eigv = ");
for(int i=0; i < n; i++){
printf("%f , ", alphai[i] / beta[i] );
}
printf("\n");
}
!A = numpy.array([[12.0, 28.0, 76.0, 220.0], [16.0, 32.0, 80.0, 224.0], [24.0, 40.0, 88.0, 232.0], [40.0, 56.0, 104.0, 248.0]], dtype='float64')
! B = numpy.array([[2.0, 4.0, 10.0, 28.0], [3.0, 5.0, 11.0, 29.0], [5.0, 7.0, 13.0, 31.0], [9.0, 11.0, 17.0, 35.0]], dtype='float64')
! D, V = scipy.linalg.eig(A, B); D
implicit none
integer, parameter :: n = 4
integer :: lda, ldb, ldvr, ldvl, lwork, info
character*1 :: jobvl, jobvr
real*8 :: alphar(n)
real*8 :: alphai(n)
real*8 :: beta(n)
real*8 :: vl(n, n)
real*8 :: vr(n, n)
real*8, allocatable :: work(:)
real*8 :: a(n, n)
real*8 :: b(n, n)
a(1, :) = (/12.0, 28.0, 76.0, 220.0/)
a(2, :) = (/16.0, 32.0, 80.0, 224.0/)
a(3, :) = (/24.0, 40.0, 88.0, 232.0/)
a(4, :) = (/40.0, 56.0, 104.0, 248.0/)
b(1, :) = (/2.0, 4.0, 10.0, 28.0/)
b(2, :) = (/3.0, 5.0, 11.0, 29.0/)
b(3, :) = (/5.0, 7.0, 13.0, 31.0/)
b(4, :) = (/9.0, 11.0, 17.0, 35.0/)
lda = n
ldb = n
ldvr = n
ldvl = n
jobvr = 'V'
jobvl = 'V'
! workspace query
allocate(work(1:8*n)) ! min value
lwork = -1
print*, 'workspace query: lwork = ', lwork
call dggev(jobvl, jobvr, n, a, lda, b, ldb, alphar, alphai, beta, vl, ldvl, vr, ldvr, work, lwork, info)
print*, 'info = ', info
lwork = int(work(1))
print*, 'opt lwork =', lwork
! do the work
deallocate(work)
allocate(work(1:lwork))
call dggev(jobvl, jobvr, n, a, lda, b, ldb, alphar, alphai, beta, vl, ldvl, vr, ldvr, work, lwork, info)
print*
print*, 'info = ', info
print*, 'alphar = ', alphar
print*, 'alphai = ', alphai
print*, 'beta = ', beta
print*
print*, 'Re(eigv) = ', alphar / beta
print*, 'Im(eigv) = ', alphai / beta
deallocate(work)
end
现在我们需要在本地编译并运行这个复制器。如果我们直接调用编译器，我们需要手动添加所需的编译和链接标志。包含路径将取决于您的本地安装，而链接标志将取决于您正在测试的库。例如，要针对 OpenBLAS 的本地构建进行测试
$ gcc ggev_repro_gh_11577.c \
-I$PWD/../flexiblas-setup/built-libs/include/ \
-L$PWD/../flexiblas-setup/built-libs/lib -lopenblas
$ ./a.out # to run the reproducer
$ gfortran ggev_repro_gh_11577.f90 \
-I/$PWD/../flexiblas-setup/built-libs/include/ \
-L$PWD/../flexiblas-setup/built-libs/lib -lopenblas
$ ./a.out # to run the reproducer
对于参考 BLAS/LAPACK，-lopenblas 应替换为 -lblas -llapack。
请注意，只有当库位于非标准位置时（例如，我们在本指南中构建的库），才需要显式路径。对于针对软件包管理器安装的库进行构建（共享库和头文件位于正常的编译器搜索路径中），例如在 /usr/lib 和 /usr/include 中，或在使用来自同一环境的编译器时在 conda 环境中，可以省略它们
$ gcc ggev_repro_gh_11577.c -lopenblas
$ ./a.out # to run the reproducer
$ gfortran ggev_repro_gh_11577.f90 -lopenblas
$ ./a.out # to run the reproducer
或者（可能是一种更稳健的方法），使用一个小的 meson.build 文件来自动执行此操作并避免手动路径
project('repro_gh_11577', 'c')
openblas_dep = dependency('openblas')
executable('repro_c',
'ggev_repro_gh_11577.c',
dependencies: openblas_dep
)
然后构建测试并运行它
$ export PKG_CONFIG_PATH=~/code/tmp/flexiblas-setup/built-libs/lib/pkgconfig/
$ meson setup build
$ ninja -C build
$ ./build/repro_c # output may vary
info = 0
Re(eigv) = 4.000000 , 8.000000 , inf , -inf ,
Im(eigv = 0.000000 , 0.000000 , -nan , -nan ,
project('repro_gh_11577', 'fortran')
openblas_dep = dependency('openblas')
executable('repro_f90',
'ggev_repro_gh_11577.f90',
dependencies: openblas_dep
)
然后构建测试并运行它
$ export PKG_CONFIG_PATH=~/code/tmp/flexiblas-setup/built-libs/lib/pkgconfig/
$ meson setup build
$ ninja -C build
$ ./build/repro_f90 # output may vary
workspace query: lwork = -1
info = 0
opt lwork = 156
info = 0
alphar = 1.0204501477442456 11.707793036240817 3.7423579363517347E-014 -1.1492523608519701E-014
alphai = 0.0000000000000000 0.0000000000000000 0.0000000000000000 0.0000000000000000
beta = 0.25511253693606051 1.4634741295300704 0.0000000000000000 0.0000000000000000
Re(eigv) = 4.0000000000000142 8.0000000000001741 Infinity -Infinity
Im(eigv) = 0.0000000000000000 0.0000000000000000 NaN NaN
警告
当您的机器上有同一 BLAS 库的多个版本/构建时，很容易在构建过程中意外地选择错误的库（请记住：-lopenblas 仅表示“链接到某个 libopenblas.so”）。如果您不确定，请对您构建的测试可执行文件使用 ldd 来检查它链接到哪个共享库。
使用 gdb 调试线性代数问题#
在调试线性代数问题时，有时需要逐步执行 Python 和 C 代码。您可以使用 pdb 进行前者，并使用 gdb 进行后者。
首先，准备一个小型的 Python 复制器，并设置一个断点。例如
$ cat chol.py
import numpy as np
from scipy import linalg
n = 40
rng = np.random.default_rng(1234)
x = rng.uniform(size=n)
a = x[:, None] @ x[None, :] + np.identity(n)
breakpoint() # note a breakpoint
linalg.cholesky(a)
然后，您需要在 gdb 下运行它，并在 LAPACK 函数处添加一个 C 级断点。这样，您的执行将停止两次：首先在 Python 断点处停止，然后在 C 断点处停止，并且您将能够逐步执行并检查 Python 和 C 变量的值。
要找出 LAPACK 名称，请读取 SciPy 函数的 Python 源代码，并在 .so 库上使用 nm 来找出确切的名称。对于上面的 Cholesky 分解，LAPACK 函数是 ?potrf，Ubuntu Linux 上的 C 名称是 dpotrf_（它可能会带有或不带有尾随下划线，大写或小写，具体取决于系统）。
这是一个 gdb 会话示例
$ gdb --args python chol.py
...
(gdb) b dpotrf_ # this adds a C breakpoint (type "y" below)
Function "dpotrf_" not defined.
Make breakpoint pending on future shared library load? (y or [n]) y
Breakpoint 1 (dpotrf_) pending.
(gdb) run # run the python script
...
> /home/br/temp/chol/chol.py(12)<module>()
-> linalg.cholesky(a) # execution stopped at the python breakpoint
(Pdb) p a.shape # ... inspect the python state here
(40, 40)
(Pdb) c # continue execution until the C breakpoint
Thread 1 "python" hit Breakpoint 1, 0x00007ffff4c48820 in dpotrf_ ()
from /home/br/miniforge/envs/scipy-dev/lib/python3.10/site-packages/numpy/core/../../../../libcblas.so.3
(gdb) s # step through the C function
Single stepping until exit from function dpotrf_,
which has no line number information.
f2py_rout__flapack_dpotrf (capi_self=<optimized out>, capi_args=<optimized out>,
capi_keywds=<optimized out>, f2py_func=0x7ffff4c48820 <dpotrf_>)
at scipy/linalg/_flapackmodule.c:63281
....
(gdb) p lda # inspect values of C variables
$1 = 40
# print out the C backtrace
(gdb) bt
#0 0x00007ffff3056b1e in f2py_rout.flapack_dpotrf ()
from /path/to/site-packages/scipy/linalg/_flapack.cpython-311-x86_64-linux-gnu.so
#1 0x0000555555734323 in _PyObject_MakeTpCall (tstate=0x555555ad0558 <_PyRuntime+166328>,
callable=<fortran at remote 0x7ffff40ffc00>, args=<optimized out>, nargs=1,
keywords=('lower', 'overwrite_a', 'clean'))
at /usr/local/src/conda/python-3.11.8/Objects/call.c:214
...
根据您的系统，您可能需要使用调试构建类型构建 SciPy，并取消设置 CFLAGS/CXXFLAGS 环境变量。有关详细信息，请参阅 NumPy 调试指南。