OPEN CASCADE Gauss Least Square
eryar@163.com
Abstract. The least square can be used to solve a set of n linear equations of m unknowns(n >= m). The OPEN CASCADE class math_GaussLeastSquare implements the least square solution of the linear equations by using Gauss LU decomposition algorithm. The paper focus on the Least Square method to solve the linear equations.
Key Words. Least Square, LU Decomposition, Linear Equations
1.Introduction
最小二乘（Least Square）問(wèn)題是一類(lèi)特殊的無(wú)約束優(yōu)化問(wèn)題，它在科學(xué)與工程計(jì)算中有十分重要的應(yīng)用。最小二乘問(wèn)題產(chǎn)生于數(shù)據(jù)擬合問(wèn)題，它是一種基于觀測(cè)數(shù)據(jù)與模型數(shù)據(jù)之間的差的平方和最小來(lái)估計(jì)模型參數(shù)的方法。它最早由德國(guó)數(shù)學(xué)家高斯Gauss于1794年，在預(yù)測(cè)行星軌道時(shí)提出，當(dāng)時(shí)高斯只有17歲，后來(lái)得到廣泛應(yīng)用。
許多工程問(wèn)題常常需要根據(jù)兩個(gè)變量的幾組實(shí)驗(yàn)數(shù)值來(lái)找出這兩個(gè)變量的函數(shù)關(guān)系的近似表達(dá)式，通常把這樣得到的函數(shù)的近似表達(dá)式稱(chēng)為經(jīng)驗(yàn)公式。經(jīng)驗(yàn)公式建立后就可以把生產(chǎn)或?qū)嶒?yàn)中所積累的某些經(jīng)驗(yàn)提高到理論上加以分析。
在幾何造型中常常需要對(duì)曲線(xiàn)和曲面進(jìn)行擬合（插值和逼近），根據(jù)一些采樣點(diǎn)來(lái)擬合曲線(xiàn)曲面。逼近比插值更為困難。在插值問(wèn)題中，只是根據(jù)采樣點(diǎn)來(lái)建立方程組，直接求解方程組即可得到結(jié)果，不需要進(jìn)行容差檢查。而在逼近問(wèn)題中，容差和采樣點(diǎn)一起作為輸入，一般預(yù)先不知道需要多少個(gè)控制點(diǎn)才能達(dá)到預(yù)期的精度，因此逼近一般都需要通過(guò)迭代來(lái)實(shí)現(xiàn)。通過(guò)最小二乘法即可實(shí)現(xiàn)達(dá)到精度要求的擬合結(jié)果，如OPEN CADCADE中的曲線(xiàn)曲面逼近就采用了最小二乘算法。
本文主要關(guān)注于最小二乘法求解線(xiàn)性方程組的原理及OPEN CASCADE中的實(shí)現(xiàn)和用法，為探索最小二乘法在OPEN CASCADE曲線(xiàn)曲面擬合方面的應(yīng)用提前做些熱身準(zhǔn)備。最小二乘問(wèn)題涉及到非線(xiàn)性最優(yōu)化的相關(guān)知識(shí)，對(duì)多元函數(shù)的微積分有些要求，可以找出原來(lái)的《高等數(shù)學(xué)》或《數(shù)學(xué)分析》的課本復(fù)習(xí)下。本文關(guān)注的應(yīng)用最小二乘法求解線(xiàn)性方程組問(wèn)題只涉及到線(xiàn)性代數(shù)或矩陣相關(guān)的知識(shí)。
2.Principle
在羅家洪、方衛(wèi)東編著的《矩陣分析引論》一書(shū)中的2.5節(jié)點(diǎn)到子空間距離與最小二乘法，用歐氏空間的概念來(lái)表達(dá)最小二乘法，并給出最小二乘解所滿(mǎn)足的代數(shù)條件的證明過(guò)程。本文摘抄主要內(nèi)容對(duì)最小二乘法求解線(xiàn)性方程組的理解。
設(shè)已給不相容實(shí)系數(shù)線(xiàn)性方程組（即無(wú)解的線(xiàn)性方程組）：
因?yàn)檫@方程組無(wú)解，設(shè)法找出一組數(shù)x1’, x2’, ..., xn’使平方偏差最小：
這組數(shù)稱(chēng)為此方程組的最小二乘解，這一方法叫做最小二乘法。經(jīng)證明，最小二乘解所滿(mǎn)足的代數(shù)方程為：
它是一個(gè)線(xiàn)性方程組，系數(shù)矩陣為ATA，常數(shù)項(xiàng)為ATB。使用上述結(jié)論來(lái)解如下線(xiàn)性方程組：
由于：
所以：
于是求得最小二乘解為：x1=17/6, x2=-13/6, x3=-4/6。
在OPEN CASCADE的數(shù)據(jù)工具集中TKMath，使用類(lèi)math_GaussLeastSquare來(lái)利用最小二乘法來(lái)對(duì)線(xiàn)性方程組進(jìn)行求解。其實(shí)現(xiàn)代碼如下所示：
math_GaussLeastSquare::math_GaussLeastSquare (const math_Matrix& A,
const Standard_Real MinPivot) :
LU(1, A.ColNumber(),
1, A.ColNumber()),
A2(1, A.ColNumber(),
1, A.RowNumber()),
Index(1, A.ColNumber()) {
A2 = A.Transposed();
LU.Multiply(A2, A);
Standard_Integer Error = LU_Decompose(LU, Index, D, MinPivot);
Done = (!Error) ? Standard_True : Standard_False;
}
void math_GaussLeastSquare::Solve(const math_Vector& B, math_Vector& X) const{
StdFail_NotDone_Raise_if(!Done, " ");
Standard_DimensionError_Raise_if((B.Length() != A2.ColNumber()) ||
(X.Length() != A2.RowNumber()), " ");
X.Multiply(A2, B);
LU_Solve(LU, Index, X);
return;
}
結(jié)合上述公式，再來(lái)理解這個(gè)代碼實(shí)現(xiàn)的思路還是很清晰的。
3.Code Example
OPEN CASCADE的TKMath工具集中提供了類(lèi)math_GaussLeastSquare實(shí)現(xiàn)了使用高斯LU分解算法求m個(gè)未知數(shù)的n個(gè)線(xiàn)性方程組的最小二乘解，其中n>=m。下面給出使用類(lèi)math_GaussLeastSquare對(duì)上述線(xiàn)性方程組進(jìn)行求解的示例程序：
/*
* Copyright (c) 2015 Shing Liu All Rights Reserved.
*
* File : main.cpp
* Author : Shing Liu(eryar@163.com)
* Date : 2015-11-25 21:00
* Version : OpenCASCADE6.9.0
*
* Description : Test Gauss Least Square method to
* solve linear equations.
*/
#define WNT
#include <math_GaussLeastSquare.hxx>
#pragma comment(lib, "TKernel.lib")
#pragma comment(lib, "TKMath.lib")
void testLeastSquare(void)
{
math_Matrix A(1, 4, 1, 3);
math_Vector B(1, 4);
math_Vector X(1, 3);
A(1,1) = 1.0; A(1,2) = 1.0; A(1,3) = 0.0; B(1) = 1.0;
A(2,1) = 1.0; A(2,2) = 0.0; A(2,3) = 1.0; B(2) = 2.0;
A(3,1) = 1.0; A(3,2) = 1.0; A(3,3) = 1.0; B(3) = 0.0;
A(4,1) = 1.0; A(4,2) = 2.0; A(4,3) = -1.0; B(4) = -1.0;
math_GaussLeastSquare aSolver(A);
aSolver.Solve(B, X);
if (aSolver.IsDone())
{
std::cout << aSolver;
std::cout << X;
}
}
int main(int argc, char* argv[])
{
testLeastSquare();
return 0;
}
程序運(yùn)行結(jié)果如下圖所示：
由上圖可知，計(jì)算結(jié)果吻合。
4.Conclusion
最小二乘法在系統(tǒng)理論中處理最小優(yōu)化問(wèn)題時(shí)有重要應(yīng)用，本文主要關(guān)注于線(xiàn)性方程組的最小二乘法求解，且對(duì)方程個(gè)數(shù)與未知數(shù)個(gè)數(shù)不要求相等。最小二乘法也是在我們學(xué)習(xí)高等數(shù)學(xué)的多元函數(shù)微分后，提出的一個(gè)實(shí)用的函數(shù)公式擬合方法。雖然本文所述的最小二乘法主要用于方程組的求解，但是OPEN CASCADE中曲線(xiàn)曲面的逼近也是采用了最小二乘法，這里最小二乘法就涉及到非線(xiàn)性最優(yōu)化的相關(guān)理論。
縱觀OPEN CASCADE的數(shù)學(xué)工具集TKMath中，大量地用到了非線(xiàn)性最優(yōu)化理論，如類(lèi)math_BFGS就實(shí)現(xiàn)了Broyden-Fletcher-Goldfarb-Shanno（BFGS），用于計(jì)算多變量函數(shù)的最小值，類(lèi)math_FRPR實(shí)現(xiàn)了Fletcher-Reeves-Polak-Ribiere算法。BFGS算法是擬牛頓方法，是解決無(wú)約束優(yōu)化問(wèn)題既快又穩(wěn)定的算法。這些最優(yōu)化算法廣泛地用于OPEN CASCADE中曲線(xiàn)曲面擬合、光順及求交等算法中。所以有必要對(duì)最優(yōu)化方法，非線(xiàn)性最優(yōu)化理論等知識(shí)進(jìn)行學(xué)習(xí)。掌握一些最優(yōu)化方法，不僅可以方便理解OPEN CASCADE中的核心關(guān)鍵算法，還可以將這些理論方法靈活應(yīng)用在自己的程序中，提高軟件質(zhì)量。由于本人能力有限，先在這兒拋磚引玉，感興趣的讀者可以結(jié)合相關(guān)書(shū)籍對(duì)非線(xiàn)性最優(yōu)化理論進(jìn)行學(xué)習(xí)，研究，應(yīng)用，創(chuàng)新。
5.References
1. 同濟(jì)大學(xué)數(shù)學(xué)教研室. 高等數(shù)學(xué). 高等教育出版社. 1996
2. 王仁宏. 李崇君. 朱春鋼. 計(jì)算幾何教程. 科學(xué)出版社. 2008
3. 羅家洪. 方衛(wèi)東. 矩陣分析引論. 華南理工大學(xué)出版社. 2006
4. 易大義. 陳道琦. 數(shù)值分析引論. 浙江大學(xué)出版社. 1998
5. 趙罡. 穆國(guó)旺. 王拉柱. 非均勻有理B樣條. 清華大學(xué)出版社. 2010
6. 王宜舉. 修乃華. 非線(xiàn)性最優(yōu)化理論與方法. 科學(xué)出版社. 2012