偵錯線性代數相關問題#
線性代數相關的錯誤報告是最難診斷和解決的問題之一。這不僅是因為線性代數在數學/演算法上可能具有挑戰性（這對於 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 這樣的發行版可能會針對 stub 函式庫（通常是 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 函式庫開發人員更容易解決問題。
要從使用 scipy 函式並以 NumPy 陣列作為輸入的 Python 重現器，轉換為 C/Fortran 重現器，有必要找到 SciPy 中採取的程式碼路徑，並確定調用了哪個確切的 BLAS 或 LAPACK 函式，以及使用什麼輸入（注意：答案可能包含在 .pyf.src f2py 簽名檔案中；查看建置目錄中生成的 _flapackmodule.c 也可能有用）。 然後可以在 C/Fortran 中通過定義一些整數/浮點變數和陣列來重現這一點（通常帶有硬編碼數字的小陣列就足夠了）。
BLAS 和 LAPACK 函式的參數列表可以在例如 Netlib LAPACK 文件 或 Reference-LAPACK/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
現在我們需要在本地編譯此重現器並執行它。 如果我們直接調用編譯器，則需要手動新增所需的編譯和連結標誌。 include 路徑將取決於您的本地安裝，而連結標誌將取決於您要測試的函式庫。 例如，要針對 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 偵錯 linalg 問題#
在偵錯 linalg 問題時，有時逐步執行 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 偵錯指南。