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

宋叢威 張曉明

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

關(guān)鍵詞：PCA; SVD; 線性方程組; 最小二乘解

文章編號：2095-2163（2019）04-0091-05 中圖分類號：TP399 文獻(xiàn)標(biāo)志碼：A

0 引 言

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

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

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

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

1 算法設(shè)計

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

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

1.1 矩陣分析

為了避免多余的矩陣運(yùn)算， 沒有按照標(biāo)準(zhǔn)的PCA流程， 而是直接對ATA進(jìn)行奇異值分解 （SVD） 或特征值分解：

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

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

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

1.2 算法

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

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

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

1.3 解的處理

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

實(shí)際上， 為了方便計算， 還可對B進(jìn)行降維，此時方程變?yōu)椋?/p>

2 數(shù)值實(shí)驗

本文算法用 Python實(shí)現(xiàn)， 運(yùn)行于MacOS上。所用數(shù)據(jù)、源代碼和運(yùn)行結(jié)果都托管在Github 上， 網(wǎng)址為https：//github.com/Freakwill/PCA-linear-equations.

2.1 主成分選擇

對矩陣A，B降維， 通過觀察累計百分比曲線， 選擇適當(dāng)?shù)闹鞒煞謹(jǐn)?shù)。如圖1所示。

2.2 誤差分析

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

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

由圖2結(jié)果可知：

（1）原方程組沒有解。圖中原方程誤差是理論上最?。ㄔ跊]有約束的情況下）， 不妨看做一個相對0點(diǎn)。如果降維方法的誤差低于這個點(diǎn)， 可認(rèn)為是系統(tǒng)誤差造成的， 不是數(shù)值分析意義上的誤差。

（2）降維方程組的誤差衰減達(dá)到預(yù)期?！柏?fù)數(shù)值置0”的處理， 對誤差沒有太大影響。當(dāng)A的主成分?jǐn)?shù)超過70時， 計算變得不穩(wěn)定， 預(yù)測誤差表現(xiàn)很夸張。

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

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

由圖4結(jié)果可知：

（1）B主成分?jǐn)?shù)對誤差影響沒有A大 （這是優(yōu)點(diǎn)）， 由于B維數(shù)較低， 且主成分比較集中， 在第一個主成分處， 誤差已降到0.3左右。

（2）用時與主成分?jǐn)?shù)近似呈弱線性相關(guān)， B的降維也相應(yīng)提高了效率。

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

可通過設(shè)定隨機(jī)種子， 產(chǎn)生固定的訓(xùn)練-測試樣本（實(shí)驗重復(fù)50次）。

選用A前30個主成分，B前4個主成分。 獲得實(shí)驗結(jié)果見表2。

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

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

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

當(dāng)主成分增加時， 負(fù)數(shù)解都是無法避免的。 數(shù)值實(shí)驗表明主成分過多還有可能出現(xiàn)異常。幾種求解方法的測試比較結(jié)果如圖5所示。

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

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

3 結(jié)束語

PCA降維在解大型線性方程組表現(xiàn)的非常出色， 隨著主成分增加， 誤差快速遞減。 這種簡單的代數(shù)學(xué)技巧， 不僅計算快， 算法設(shè)計簡單， 由于矩陣分解為對角矩陣和正交矩陣的乘積， 之后也無需解任何線性方程組?？梢愿鶕?jù)需要任意壓縮原方程。

本文提供的方法可以應(yīng)用于工業(yè)生產(chǎn)中。但在實(shí)驗中發(fā)現(xiàn)過多的主成分可能使得算法不穩(wěn)定。今后的工作重點(diǎn)：

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

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

（3）應(yīng)利用系數(shù)矩陣的一些特點(diǎn)， 如稀疏性， 設(shè)計更有針對性的算法。

參考文獻(xiàn)

[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] 焦李成，等. 稀疏學(xué)習(xí)、分類與識別[M]. 北京：科學(xué)出版社， 2017.

[6] 解鋒昌， 韋博成， 林金官. 零過多數(shù)據(jù)的統(tǒng)計分析及其應(yīng)用[M]. 北京：科學(xué)出版社， 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/.