宋叢威 張曉明

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

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

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

0 引 言

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

目前，印染行業(yè)急需解決的難題是在染色過(guò)程中， 因?yàn)槭孪炔磺宄玖瞎庾V的精確數(shù)值， 只能先從多個(gè)布匹的染色流程中， 獲取染色光譜數(shù)據(jù)， 再估計(jì)使用染料的光譜數(shù)據(jù)， 最后用于新的布匹染色。設(shè)建立一個(gè)線性模型：

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

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

1 算法設(shè)計(jì)

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

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

1.1 矩陣分析

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

由式 （6） 可知， 對(duì)誤差來(lái)說(shuō)， 選擇哪些主成分似乎并不重要， 相對(duì)誤差大致可以寫(xiě)成：

若對(duì)B降維，有如下公式：

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

1.2 算法

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

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

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

1.3 解的處理

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

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

2 數(shù)值實(shí)驗(yàn)

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

2.1 主成分選擇

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

2.2 誤差分析

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

20%的樣本會(huì)被事先抽取出來(lái)， 測(cè)試誤差通常比訓(xùn)練誤差更有說(shuō)服力。降維是對(duì)訓(xùn)練樣本進(jìn)行的， 而不是對(duì)測(cè)試樣本進(jìn)行。但對(duì)A的降維也可以應(yīng)用于訓(xùn)練樣本， 因?yàn)槠涫且阎斎胫担?而測(cè)試樣本中的B作為目標(biāo)值， 不參與降維。如圖2所示。

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

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

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

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

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

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

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

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

總之， 對(duì)B降維是完全合理的。

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

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

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

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

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

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

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

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

3 結(jié)束語(yǔ)

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

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

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

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

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

參考文獻(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í)、分類(lèi)與識(shí)別[M]. 北京：科學(xué)出版社， 2017.

[6] 解鋒昌， 韋博成， 林金官. 零過(guò)多數(shù)據(jù)的統(tǒng)計(jì)分析及其應(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/.