【Nature Methods 】SciPy 1.0:Python科學計算的基本算法

       SciPy是Python程式語言的一個數學程序庫，它為建模和解決科學問題提供了構建塊。SciPy包括優化、積分、插值、特徵值問題、代數方程、微分方程和許多其他類問題的算法；它還提供特殊的數據結構，如稀疏矩陣和k-維樹。SciPy是建立在NUMPY之上的，它提供了array數據結構和相關的快速的數學程序，而SciPy本身就是構建更高級別的科學庫的基礎，包括scikit-learn和scikit-image。在這篇文章中，我們概述了SciPy 1.0的功能和開發實踐，並重點介紹了最近的一些技術進展。

    截至2019年2月，SciPy庫包含將近600,000行開原始碼，這些代碼按圖1匯總的16個子包進行組織。開發團隊和社區目前主要在GitHub上進行交互和操作，這是一個在線版本控制和任務管理平臺。超過110,000個GitHub存儲庫和6,500個軟體包依賴SciPy。

圖1：SciPy包組織

    一、架構和實現

    1.1 項目範圍

    SciPy提供了用於科學計算的基本算法。在發展相對較緩慢的領域，SciPy旨在提供完整的覆蓋範圍。在其他領域，SciPy的目的是提供基本的構建塊，同時與該領域的其他包進行良好的交互。例如，SciPy提供了概率分布，假設檢驗，頻率統計，相關函數等，而Statsmodels提供了更高級的統計估計量和推斷方法，scikit- learn3涵蓋了機器學習，而PyMC3，emcee和PyStan  （http://mc-stan.org） 涵蓋了貝葉斯統計和概率建模。scikit-image4提供了比SciPy的ndimage模塊更加優秀的的圖像處理功能，SymPy提供了用於符號計算的Python交互界面，而sparse. csgraph和spatial與諸如NetworkX之類的專用庫相比提供了與圖形和網絡一起使用的基本工具。

    我們使用以下標準來確定是否要在SciPy中加入新特徵：

   

     •該算法與多個科學領域相關。

   

     •該算法非常重要。例如，它足夠經典，可以包含在教科書中，或者基於經過大量同行引用的同行評審文章。

    

    在軟體系統和體系結構方面，SciPy的範圍與NumPy的範圍相匹配。

 

    1.2 語言選擇

    根據linguist庫 （https://github.com /github/linguist） 的分析，SciPy大約是Python佔50％，Fortran佔25％，C佔20％，Cython佔3％和C ++佔2％，以及少量的TeX， Matlab，shell 腳本和Make。

    儘管Fortran已有很長的歷史，但它仍然是一種高性能的科學程式語言。因此，我們包裝了以下出色的且經過實測的Fortran庫，以在為Python帶來便利的同時提高性能：QUADPACK和ODEPACK用於數值積分和初值問題的解決；FITPACK，ODRPACK和MINPACK用於曲線擬合和最小二乘最小化；FFTPACK用於執行傅立葉變換；ARPACK用於解決特徵值問題；用於計算貝塞爾函數的ALGORITHM 644；和CDFLIB用於評估累積密度函數。

    在SciPy中排名前三位的語言還有一個是C語言，我們在SciPy中包裝的C庫包括：用於優化的trlib，用於解決稀疏線性系統的SuperLU，用於計算幾何圖形的Qhull和用於特殊函數的Cephes （http://www.netlib.org /cephes/）。

    Cython被描述為一種混合語言，它混合了Python的最佳部分和低級C/C++範例。我們經常使用Cython作為C/C++編寫的成熟的低級科學計算庫與SciPy提供的Python接口之間的粘合劑。我們還使用Cython在Python代碼中實現性能增強，特別是對於大量使用的內部循環受益於具有靜態類型的編譯代碼的情況。

    為了實現新功能，Python仍然是首選的語言。如果Python性能存在問題，那麼我們更喜歡使用Cython，然後使用C，C++或Fortran（按此順序）。其主要考慮因素是可維護性：Cython具有最高的抽象級別，大多數Python開發人員都會理解它。C也是很常用的，比起C++尤其是Fortran，當前的核心開發團隊更容易管理。

 

    1.3 API和ABI的演變

    SciPy的API由大約1500個函數和類組成。我們優化API的策略是添加新功能，而只有在優點超過用戶的損失以及向這些用戶提供明確的捨棄警告至少一年之後才能刪除或更改現有功能。一般而言，我們鼓勵修改從而提高庫的API的清晰度，但不建議破壞向後兼容性，因為SciPy位置接近科學Python計算堆棧的基礎。

    

    除了Python API之外，SciPy還有C和Cython接口。因此，我們還必須考慮應用程式二進位接口（ABI）。這個ABI已經穩定了很長時間，我們的目標是以向後兼容的方式發展它。

 

    二、關鍵技術改進

    在這裡，我們描述了過去三年中取得的關鍵技術改進。

    2.1 數據結構

    稀疏矩陣。scipy. sparse提供了七個稀疏矩陣數據結構，也稱為稀疏格式。最重要的是行和列壓縮格式（分別為CSR和CSC）。它們提供快速的軸索引和快速的矩陣向量乘法，在整個SciPy和相關程序包中大量使用。

    在過去三年中，我們對稀疏矩陣處理內部進行了重寫，提高了性能。現在，CSC和CSR矩陣的迭代和切片速度最高可提高35％，並且COO/DIA轉換為CSR/ CSC格式的速度有所提高。SuperLU已更新至版本5.2.1，增強了部分稀疏矩陣功能的低級實現。

    從新功能的角度來看，scipy.sparse矩陣和線性運算操作現在支持Python矩陣乘法（@）運算符。我們添加了scipy. sparse. norm和scipy. sparse. random分別用於計算稀疏矩陣範數和從任意分布中得出隨機變量。此外，我們盡力使scipy. sparse API與NumPy API保持一致。

 

    cKDTree。scipy. spatial. ckdtree模塊實現了組織k維空間中的點的空間分區數據結構，並使用模板化類在C++中進行了重寫。增加了對周期性邊界條件的支持，該條件通常用於物理過程的仿真中。

    2013年，cKDTree.query的k最近鄰搜索的時間複雜度大約為loglinear，與其正式描述一致。從那時起，我們增強了cKDTree,通過在C ++中重新實現該查詢，消除內存洩漏並允許釋放全局解釋器鎖（GIL），以便可以使用多個線程來進行查詢。這通常可以在保留漸近複雜性的同時提高任何給定問題的性能。

    在2015年，SciPy添加了sparse_ distance _matrix例程，通過忽略所有超過用戶提供的值的距離來生成KDTree對象之間的近似稀疏距離矩陣。該例程不限於常規的L2（歐幾裡得）範數，而是支持1到無窮大之間的任何Minkowski p-norm。默認情況下，返回的數據結構是基於字典（DOK）的稀疏矩陣，這對於矩陣構造非常有效。這種用於稀疏矩陣彙編的哈希方法可以比使用CSR 格式構造快7倍，並且C++級別的稀疏矩陣構造釋放了Python GIL以提高性能。

    2015年，對cKDTree算法進行了增強，以支持權重，這在許多科學應用中都是必不可少的，例如，計算星系的相關函數。

 

    2.2 對已編譯代碼的統一綁定

    LowLevelCallable。從SciPy 0.19版開始，用戶可以將低級函數包裝在scipy. LowLevelCallable對象中，從而減少了直接從Python調用已編譯的C函數（例如使用Numba或Cython生成的函數）的開銷。支持的低級功能包括PyCapsule對象，ctypes函數指針和cffi函數指針。此外，可以從Cython模塊使用scipy. LowLevelCallable. from _ cython自動生成低級回調函數。

 

    2.3 BLAS，LAPACK和special的Cython綁定

    SciPy現在包括許多BLAS和LAPACK程序以及scipy. special子包中提供的特殊函數的Cython包裝器。這些用於Cython的低級接口也可以在SciPy代碼庫之外使用，以訪問包裝的庫中的函數，同時避免Python函數調用的開銷。對於許多用例，這可以使性能提高一兩個數量級。開發人員也可以使用低級Cython接口，而無需連結到所包裝的庫。這使其他擴展功能避免了查找和使用正確庫的複雜性。

    包裝以Fortran編寫的庫時，避免這種複雜性尤為重要。這些低級包裝器不僅可以在沒有Fortran編譯器的情況下使用，而且還可以在無需處理所有不同的Fortran編譯器ABI和名稱修改方案的情況下使用。這些低級Cython包裝器大多數都是自動生成的，以保證正確性和易於維護性。

    由於SciPy可以使用LAPACK 3.4.0或更高版本構建，因此僅為一致的程序提供Cython包裝器所有支持的LAPACK版本之間的接口。各種現有BLAS庫提供的標準BLAS接口當前未更改，因此SciPy提供的包裝程序通常不需要更改。scipy. special的Cython包裝器遵循對該子包的接口的更改進行相應更改。

 

    2.4 數值優化

    scipy. optimize子程序包提供了一些用於求根和優化問題函數。在這裡，我們重點介紹SciPy1.0的最新功能。

    線性優化。與SciPy 1.0一起發布了一種新的針對連續線性編程問題的內點優化器linprog，method = 』interior-point』,實現了商用求解器MOSEK的核心算法，它可以解決所有90多個經過測試的NETLIB LP基準測試問題。

    預求解程序解決了一些瑣碎的問題，在形式上簡化了問題，例如束縛收緊和移除固定變量，並且自動選擇了一些用於消除冗餘等式約束的程序，以減少由奇異矩陣引起的數值計算困難的可能性。儘管主要的求解器是由Python實現，但由於端到端稀疏矩陣的支持和SciPy編譯的線性系統求解器的大量使用，可提供足夠的速度來解決有成千上萬的變量和約束的問題。

    非線性優化：局部最小化。minimize功能提供了一個統一的接口，用於查找非線性優化問題的局部最小值。在最新版本的SciPy中，在minimize中添加了四種用於無約束優化的新方法：dogleg，trust-ncg，trust-exact和trust-krylov。所有方法都是信賴域方法，它們基於一階和二階導數信息建立目標函數的局部模型，逼近局部「信賴域」內的最佳點，並進行迭代，直到達到原始目標函數的局部最小值為止，但是每種方法都有其獨特的特徵，使其適合某些類型的問題。例如，trust-exact通過幾乎精確地解決信賴區域子問題來實現快速收斂，但是它要求存儲第二階導數Hessian矩陣並在每次迭代中將其分解，這可能會影響大規模問題（≥1,000個變量）的求解。相反，trust-ncg和trust-krylov非常適合大規模優化問題，因為它們不需要顯式存儲和分解Hessian矩陣，而可以以更快的近似方式使用二階導數信息。我們在表1中詳細比較了所有最小化方法的特性，這些特性說明了SciPy涵蓋數值方法時所要達到的完整性水平。

 

表1：minimize的優化算法

    非線性優化：全局最小化。由於最小化可能返回任何局部最小值，某些問題需要使用全局最小值。新的scipy. optimize. dif - ferentialevolution函數是一種隨機全局優化器，它通過演化大量候選解而起作用。

    在每個迭代中，通過組合現有組中的候選來生成試驗候選組。如果候選試驗代表一種優化，那麼將更新組。最近，SciPy基準測試套件獲得了196個全局優化問題的綜合集合，這些問題可以隨著時間的推移跟蹤現有求解器的性能，並評估新求解器的性能是否值得將其包含在軟體包中。

 

    2.5 統計分布

    scipy. stats子程序包包含100多個概率分布：96個連續分布和13個離散單變量分布，以及10個多元分布。

    該實現依賴於一個一致的框架，該框架提供了對隨機變量進行採樣，評估累積分布函數（CDF）和概率密度函數（PDF）以及為每個分布擬合參數的方法。通常，這些方法依賴於每種分布的特定實現，例如CDF的閉式表達式或採樣算法（如果可用）。否則，將基於通用代碼使用默認方法，例如，對PDF進行數值積分以獲得CDF。最近添加到scipy. stats中的主要分布包括在scipy. stats中基於直方圖的分布,rv _ histogram和scipy中的多項式分布stats. multinomial（用於自然語言處理）

    2.6 多項式插值器

    從歷史上看，SciPy嚴重依賴於P.  Dierckx的FITPACK Fortran庫進行數據的單變量插值和逼近，但是SciPy和FITPACK之間交互的原始單片設計和API限制了用戶和開發人員。

    實現多項式插值器的一種新的模塊化設計的過程遍布多個版本。這項工作的目標是要有一組代表分段多項式的基本對象，以實現用於構建各種插值器的算法的集合，並為用戶提供用於構建其他插值器的構件。

    在新設計的最低層是代表單變量分段多項式的類：PPoly （SciPy 0.13） ，BPoly （SciPy 0.13） 和BSpline （SciPy0.19） ，它們允許進行有效的矢量化評估，微分，積分和求根。PPoly以冪為單位，以斷點和每個區間的係數表示分段多項式。BPoly是類似的，並且在Bernstein基礎上表示分段多項式（例如，構造Bézier曲線）。BSpline表示樣條曲線，即B樣條基礎元素的線性組合。

    在下一層中，這些多項式類用於構造幾種常見的數據插值方式：CubicSpline （SciPy 0.18） 構造一個二次可微的分段三次函數，Akima1DInterpolator和PCHIPInterpolator用於構造C連續單調形狀保全插值器。

 

    2.7 測試和基準套件

    測試套件。對於SciPy的每個組件，我們編寫了多個小型可執行測試，以驗證其預期的行為。根據持續集成的實踐，在運行測試套件之前，所有對SciPy的建議貢獻都將與庫的主分支暫時集成在一起，並且在永久合併貢獻之前必須通過所有測試。持續監視單元測試所涵蓋的SciPy中的代碼行數是我們可以確定更改和正確實現新功能的一種方法。

    SciPy測試套件由一個連續集成矩陣編排而成，該矩陣包括分別由Travis CI和AppVeyor管理的POSIX和Windows （32/64位） 平臺。我們的測試涵蓋Python版本2.7、3.4、3.5、3.6，並使用pyflakes和pycodestyle進行代碼替換。測試套件中有13,000多個單元測試，這些單元測試是為與pytest （https:// docs.pytest.org /en /latest） 框架一起使用而編寫的。在圖2中，我們顯示了使用基於Docker的方法 （https:// github.com/ tylerjereddy/ scipy-cov-track） 生成的歷史測試覆蓋率數據。根據pytest-cov （https:// pypi.org/ project/ pytest-cov/） ，對於Python代碼，SciPy 1.0發行點的測試覆蓋率為87％。CircleCI服務會自動生成和發布代碼文檔，以便於評估文檔更改/完整性。

 

圖2：Python和SciPy中的編譯代碼量

    基準套件。除了確保單元測試通過外，重要的是要確認SciPy代碼庫的性能會隨著時間的推移而提高。自2015年2月以來，SciPy的性能已通過Airspeed Velocity （asv https:// github.com/ airspeed-velocity/ asv） 進行監控。SciPy的run.py腳本可以方便地包裝asv功能，以便可以通過單個控制臺命令隨時間生成基準測試結果。例如，在圖3中，展示了scipy. spatial的改進。cKDTree. query大約有九年的項目歷史。在基準測試中使用的樹是在沒有應用環形拓撲的情況下生成的 （boxsize = None） ，並且使用Airspeed Velocity 0.4,Python 2.7，NumPy 1.8.2和Cython 0.27.3、0.21.、0.18版本進行了測試，向後兼容。將cKDTree完全進行Cython化後，再用C++對其進行重寫時，實現了顯著的性能改進。

圖3：從cKDTree的引入到SciPy 1.0發行的基準測試結果

    三、討論

    SciPy具有強大的開發者社區和龐大的用戶群。儘管如此，SciPy一直在努力改進。許多開源項目面臨的問題是吸引和留住開發人員。因為代碼交接過多會導致記憶喪失，從而導致過去的錯誤重複出現。幸運的是，SciPy項目繼續吸引著熱情而有才幹的新開發人員，同時保持了規模雖小但專職的前期貢獻者的參與。

    要解決的最後一個重要挑戰是如何容納GPU和分布式計算，而又不會破壞我們傳統的架構。儘管我們仍不清楚如何在整個庫中採用這些新興技術的確切方法，但是我們現在有一個具體的子包實現，該子包允許實驗性地使用多個後端，這將在以後的報告中詳細描述。

    本文來源，可見文末的 「閱讀原文」 。