基于PCA的解大型超定線性方程組快速算法及應用

宋叢威 張曉明

摘 要：本文用主成分分析 （PCA） 降維技術解大型超定線性方程組AX=B的最小二乘解。工業上遇到的大型方程組的求解要消耗大量的時間和內存， 系數矩陣通常是病態的， 會帶來不可忽視的誤差。本文設計基于PCA的算法， 并從理論上說明其可行性。實驗證明該方法有效， 不僅測試誤差極小， 接近原方程的誤差， 而且計算時間顯著減少。

關鍵詞：PCA; SVD; 線性方程組; 最小二乘解

文章編號：2095-2163（2019）04-0091-05 中圖分類號：TP399 文獻標志碼：A

0 引 言

很多工業問題最終會轉化成解線性方程組問題。但是，通常工業級的方程組規模巨大， 直接求解會消耗大量的時間和內存。本文提出一種基于主成分分析 （PCA） [1-5]的快速降維算法求解大型線性方程組。

目前，印染行業急需解決的難題是在染色過程中， 因為事先不清楚染料光譜的精確數值， 只能先從多個布匹的染色流程中， 獲取染色光譜數據， 再估計使用染料的光譜數據， 最后用于新的布匹染色。設建立一個線性模型：

N×n。問題是：使用了多種染料， 獲取了大量光譜數據。方程規模非常巨大， 常規解法可能會消耗大量時間和內存。因此用機器學習的降維方法來縮小矩陣大小， 然后再解方程。

本文采用了5 000多組數據，（而這只是工業數據中極小的一部分）， 目的在于驗證本文方法的有效性。

1 算法設計

解線性方程組AX=B， 其中矩陣A很大， 可能是病態的， 且一般是列滿秩的 （A的行數遠大于列數）， 因此可用ATA的條件數衡量方程病態程度。B有時也不只有一列。

由于各種原因， 原方程可能沒有解， 只能尋求最小二乘解， 即解優化問題。

1.1 矩陣分析

為了避免多余的矩陣運算， 沒有按照標準的PCA流程， 而是直接對ATA進行奇異值分解 （SVD） 或特征值分解：

由式 （6） 可知， 對誤差來說， 選擇哪些主成分似乎并不重要， 相對誤差大致可以寫成：

若對B降維，有如下公式：

其中，B1是對B的PCA重構。觀察式 （7）， 顯然應該選取對應特征值絕對值最大的主成分。

1.2 算法

根據上述分析設計如下算法（見表1），將一個大型方程組分解成2個小型方程組。

若使用NMF或其它分解技術， 分解A=WH， 則不能避免做最小二乘法。

理論上， 最小二乘解是方程誤差最小的解， 但降維之后， 這個解就不再是原方程的最優解了。

1.3 解的處理

設最后得到的[AKX^]=V1Y1，是原方程AX=B的降維（最小二乘）解。工業中有時候默認所有量都是非負的， 但是降維解可能含有負數。為了滿足非負性， 可將其修正為：

實際上， 為了方便計算， 還可對B進行降維，此時方程變為：

2 數值實驗

本文算法用 Python實現， 運行于MacOS上。所用數據、源代碼和運行結果都托管在Github 上， 網址為https：//github.com/Freakwill/PCA-linear-equations.

2.1 主成分選擇

對矩陣A，B降維， 通過觀察累計百分比曲線， 選擇適當的主成分數。如圖1所示。

2.2 誤差分析

誤差采用向量的2范數， 對矩陣來說就是Frobenius 范數。本文用相對誤差公式衡量算法逼近能力。

20%的樣本會被事先抽取出來， 測試誤差通常比訓練誤差更有說服力。降維是對訓練樣本進行的， 而不是對測試樣本進行。但對A的降維也可以應用于訓練樣本， 因為其是已知輸入值， 而測試樣本中的B作為目標值， 不參與降維。如圖2所示。

由圖2結果可知：

（1）原方程組沒有解。圖中原方程誤差是理論上最小（在沒有約束的情況下）， 不妨看做一個相對0點。如果降維方法的誤差低于這個點， 可認為是系統誤差造成的， 不是數值分析意義上的誤差。

（2）降維方程組的誤差衰減達到預期。“負數值置0”的處理， 對誤差沒有太大影響。當A的主成分數超過70時， 計算變得不穩定， 預測誤差表現很夸張。

（3）從誤差曲線不難看出， 解方程組時并不需要用到所有樣本，樣本數量對本算法沒有太大影響。實際應用中， 隨機挑選充分數量的樣本即可。

圖3是對B取前4個主成分得到的誤差圖， 和之前的實驗相比，并沒有顯著影響誤差。時間和主成分數基本上呈線性關系， 符合預期。為了獲得較為準確的運行時間， 本次實驗對每一個主成分數重復50次， 無論時間還是誤差都取均值。

由圖4結果可知：

（1）B主成分數對誤差影響沒有A大 （這是優點）， 由于B維數較低， 且主成分比較集中， 在第一個主成分處， 誤差已降到0.3左右。

（2）用時與主成分數近似呈弱線性相關， B的降維也相應提高了效率。

總之， 對B降維是完全合理的。

可通過設定隨機種子， 產生固定的訓練-測試樣本（實驗重復50次）。

選用A前30個主成分，B前4個主成分。 獲得實驗結果見表2。

2.3 其它降維策略與擬合方法

NMF降維的效果和PCA相近， 但矩陣分解后需解較為復雜的方程。PCA的優勢在于能分解出正交矩陣， 可以設計出更快捷的算法， 計算用時短。中心化PCA在主成分數較少時表現較好， 這是因為中心化PCA采用的是仿射變換， 比單純的線性變換多了常數項， 但隨著主成分數增加， 并沒有表現出明顯優勢。目前中心化PCA只是簡單重構A， 即：

其中，C1VT1是中心化后的分解， 即通常意義上的PCA， M是均值矩陣，其無法簡化方程， 計算時間維持在某個常數。若能設計出簡化方程的算法， 那么中心化PCA或許是不錯的選擇。

當主成分增加時， 負數解都是無法避免的。 數值實驗表明主成分過多還有可能出現異常。幾種求解方法的測試比較結果如圖5所示。

顯然，完全可以用其它方法求解降維后的方程組， 尤其是當人們只關心預測， 而不在乎X時。如果只為了建立A與B的聯系， 完全可以采用一些非線性方法。表3中給出一些線性模型測試結果，用到了一些較復雜的計算策略， 含非線性成分， 誤差無顯著降低， 運行時間則較長。這些模型均由Python 第三方庫scikit-learn實現[7-8]。

最后， 給出一般的計算框架， 可為解線性方程組開發新的算法：

3 結束語

PCA降維在解大型線性方程組表現的非常出色， 隨著主成分增加， 誤差快速遞減。 這種簡單的代數學技巧， 不僅計算快， 算法設計簡單， 由于矩陣分解為對角矩陣和正交矩陣的乘積， 之后也無需解任何線性方程組。可以根據需要任意壓縮原方程。

本文提供的方法可以應用于工業生產中。但在實驗中發現過多的主成分可能使得算法不穩定。今后的工作重點：

（1）當變量有非負性或其它約束條件時， 本文還沒有給出更合理的處理辦法。

（2）如何實現有效的基于中心化PCA的算法。本文提到的中心化PCA沒有真正起到壓縮方程組的作用。

（3）應利用系數矩陣的一些特點， 如稀疏性， 設計更有針對性的算法。

參考文獻

[1]HASTIE T， TIBSHIRANI R， FRIEDMAN J H. The elements of statistical learning：Data mining， inference， and prediction [M]. New York ：Springer-Verlag， 2001.

[2] HOTELLING H. Analysis of a complex of statistical variables into principal components[J]. Journal of Educational Psychology， 1993，24（6）：417-441.

[3] TURK M A， PENTLAND A P. Face recognition using eigenfaces[C]//Proceedings 1991 IEEE Computer Society Conference on ComputerVision and Pattern Recognition. Maui， HI， USA：IEEE，1991：586-591.

[4] SHALEV-SHWARTZ S， BEN-DAVID S. Understanding machine learning：from theory to algorithm[M]. New York， NY， USA ：Cambridge University Press， 2014.

[5] 焦李成，等. 稀疏學習、分類與識別[M]. 北京：科學出版社， 2017.

[6] 解鋒昌， 韋博成， 林金官. 零過多數據的統計分析及其應用[M]. 北京：科學出版社， 2013.

[7] HACKELING G. Mastering machine learning with scikit-learn[M]. 2nd ed. Birmingham ：Packt Publishing， 2017.

[8] COURNAPEAU D. Scikit-learn[EB/OL]. [2019]. https：//scikit-learn.org/stable/.