線性模型（二）之多項式擬合
創(chuàng)新互聯(lián)建站-專業(yè)網(wǎng)站定制、快速模板網(wǎng)站建設(shè)、高性價比海林網(wǎng)站開發(fā)、企業(yè)建站全套包干低至880元,成熟完善的模板庫,直接使用。一站式海林網(wǎng)站制作公司更省心,省錢,快速模板網(wǎng)站建設(shè)找我們，業(yè)務(wù)覆蓋海林地區(qū)。費用合理售后完善，10余年實體公司更值得信賴。
1. 多項式擬合問題
??多項式擬合（polynominal curve fitting）是一種線性模型，模型和擬合參數(shù)的關(guān)系是線性的。多項式擬合的輸入是一維的，即x=xx=x，這是多項式擬合和線性回歸問題的主要區(qū)別之一。
??多項式擬合的目標(biāo)是構(gòu)造輸入xx的MM階多項式函數(shù)，使得該多項式能夠近似表示輸入xx和輸出yy的關(guān)系，雖然實際上xx和yy的關(guān)系并不一定是多項式，但使用足夠多的階數(shù)，總是可以逼近表示輸入xx和輸出yy的關(guān)系的。
??多項式擬合問題的輸入可以表示如下：
D={(x1,y1),(x2,y2),...,(xi,yi),...,(xN,yN)}xi∈Ryi∈R
D={(x1,y1),(x2,y2),...,(xi,yi),...,(xN,yN)}xi∈Ryi∈R
??目標(biāo)輸出是得到一個多項式函數(shù)：
f(x)=w1x1+w2x2+wixi+...+wMxM+b=(∑i=1Mwixi)+b
f(x)=w1x1+w2x2+wixi+...+wMxM+b=(∑i=1Mwixi)+b
其中MM表示最高階數(shù)為MM。
??可見在線性擬合的模型中，共包括了(M+1)(M+1)個參數(shù)，而該模型雖然不是輸入xx的線性函數(shù)，但卻是(M+1)(M+1)個擬合參數(shù)的線性函數(shù)，所以稱多項式擬合為線性模型。對于多項式擬合問題，其實就是要確定這(M+1)(M+1)個參數(shù)，這里先假設(shè)階數(shù)MM是固定的（MM是一個超參數(shù)，可以用驗證集來確定MM最優(yōu)的值，詳細(xì)的關(guān)于MM值確定的問題，后面再討論），重點就在于如何求出這(M+1)(M+1)個參數(shù)的值。
2.優(yōu)化目標(biāo)
??多項式擬合是利用多項式函數(shù)逼近輸入xx和輸出yy的函數(shù)關(guān)系，通過什么指標(biāo)來衡量某個多項式函數(shù)的逼近程度呢？（其實這就是誤差/損失函數(shù)）。擬合/回歸問題常用的評價指標(biāo)是均方誤差（在機器學(xué)習(xí)中的模型評估與度量博客中，我進行了介紹）。多項式擬合問題也同樣采用該評價指標(biāo)，以均方誤差作為誤差/損失函數(shù)，誤差函數(shù)越小，模型越好。
E(w,b)=1N∑i=1N[f(xi)?yi]2
E(w,b)=1N∑i=1N[f(xi)?yi]2
??系數(shù)1N1N是一常數(shù)，對優(yōu)化結(jié)果無影響，可以去除，即將均方誤差替換為平方誤差：
E(w,b)=∑i=1N[f(xi)?yi]2
E(w,b)=∑i=1N[f(xi)?yi]2
?? 到這里，就成功把多項式擬合問題變成了最優(yōu)化問題，優(yōu)化問題可表示為：
argminw,bE(w,b)
arg?minw,b?E(w,b)
即需要求得參數(shù){w1,...,wM,b}{w1,...,wM,b}的值，使得E(w,b)E(w,b)最小化。那么如何對該最優(yōu)化問題求解呢？
3. 優(yōu)化問題求解
3.1 求偏導(dǎo)，聯(lián)立方程求解
?? 直觀的想法是，直接對所有參數(shù)求偏導(dǎo)，令偏導(dǎo)為0，再聯(lián)立這M+1M+1個方程求解（因為共有M+1M+1個參數(shù)，故求偏導(dǎo)后也是得到M+1M+1個方程）。
E(w,b)=∑i=1N[f(xi)?yi]2=∑i=1N[(w1x1i+w2x2i+wixji+...+wMxMi+b)?yi]2
E(w,b)=∑i=1N[f(xi)?yi]2=∑i=1N[(w1xi1+w2xi2+wixij+...+wMxiM+b)?yi]2
利用E(w,b)E(w,b)對各個參數(shù)求偏導(dǎo)，如下：
?E(w,b)?wj?E(w,b)?b=2∑i=1N[(w1x1i+w2x2i+wixji+...+wMxMi+b)?yi]xji=2∑i=1N[(w1x1i+w2x2i+wixji+...+wMxMi+b)?yi]
?E(w,b)?wj=2∑i=1N[(w1xi1+w2xi2+wixij+...+wMxiM+b)?yi]xij?E(w,b)?b=2∑i=1N[(w1xi1+w2xi2+wixij+...+wMxiM+b)?yi]
求導(dǎo)之后，將各個點(xi,yi)(xi,yi)的值帶入偏導(dǎo)公式，聯(lián)立方程求解即可。
??針對該解法，可以舉個例子詳細(xì)說明，比如有兩個點(2,3),(5,8)(2,3),(5,8),需要利用二階多項式f(x)=w1x+w2x2+bf(x)=w1x+w2x2+b擬合。求解過程如下：
該二階多項式對參數(shù)求偏導(dǎo)得到
?E(w,b)?wj?E(w,b)?b=2∑i=12[(w1x1i+w2x2i+b)?yi]xji=[(w1x1+w2x21+b)?y1]xj1+[(w1x2+w2x22+b)?y2]xj2=2∑i=12[(w1x1i+w2x2i+b)?yi]=[(w1x1+w2x21+b)?y1]+[(w1x2+w2x22+b)?y2]
?E(w,b)?wj=2∑i=12[(w1xi1+w2xi2+b)?yi]xij=[(w1x1+w2x12+b)?y1]x1j+[(w1x2+w2x22+b)?y2]x2j?E(w,b)?b=2∑i=12[(w1xi1+w2xi2+b)?yi]=[(w1x1+w2x12+b)?y1]+[(w1x2+w2x22+b)?y2]
將點(2,3),(5,8)(2,3),(5,8)帶入方程，可以得到3個方程，
2b+7w1+29w2=117b+29w1+133w2=4629b+133w1+641w2=212
2b+7w1+29w2=117b+29w1+133w2=4629b+133w1+641w2=212
聯(lián)立這三個方程求解，發(fā)現(xiàn)有無窮多的解，只能得到3w1+21w2=53w1+21w2=5，這三個方程是線性相關(guān)的，故沒有唯一解。
??該方法通過求偏導(dǎo)，再聯(lián)立方程求解，比較復(fù)雜，看著也很不美觀。那么有沒有更加方便的方法呢？
3.2 最小二乘法
?? 其實求解該最優(yōu)化問題（平方和的最小值）一般會采用最小二乘法（其實最小二乘法和求偏導(dǎo)再聯(lián)立方程求解的方法無本質(zhì)區(qū)別，求偏導(dǎo)也是最小二乘法，只是這里介紹最小二乘的矩陣形式而已）。最小二乘法（least squares），從英文名非常容易想到，該方法就是求解平方和的最小值的方法。
??可以將誤差函數(shù)以矩陣的表示(NN個點，最高MM階)為：
∥Xw?y∥2
‖Xw?y‖2
其中，把偏置bb融合到了參數(shù)ww中，
w={b,w1,w2,...,wM}
w={b,w1,w2,...,wM}
XX則表示輸入矩陣，
??????11...1x1x2...xNx21x22...x2N............xM1xM2...xMN??????
[1x1x12...x1M1x2x22...x2M...............1xNxN2...xNM]
yy則表示標(biāo)注向量，
y={y1,y2,...,yN}T
y={y1,y2,...,yN}T
因此，最優(yōu)化問題可以重新表示為
minw∥Xw?y∥2
minw‖Xw?y‖2
對其求導(dǎo)，
?∥Xw?y∥2?w=?(Xw?y)T(Xw?y)?w=?(wTXT?yT)(Xw?y)?w=?(wTXTXw?yTXw?wTXTy+yTy)?w
?‖Xw?y‖2?w=?(Xw?y)T(Xw?y)?w=?(wTXT?yT)(Xw?y)?w=?(wTXTXw?yTXw?wTXTy+yTy)?w
在繼續(xù)對其求導(dǎo)之前，需要先補充一些矩陣求導(dǎo)的先驗知識（常見的一些矩陣求導(dǎo)公式可以參見轉(zhuǎn)載的博客），如下：
?xTa?x=a?ax?x=aT?xTA?x=Ax+ATx
?xTa?x=a?ax?x=aT?xTA?x=Ax+ATx
根據(jù)上面的矩陣求導(dǎo)規(guī)則，繼續(xù)進行損失函數(shù)的求導(dǎo)
?∥Xw?y∥2?w=?(wTXTXw?yTXw?wTXTy+yTy)?w=XTXw+(XTX)Tw?(yTX)T?XTy=2XTXw?2XTy
?‖Xw?y‖2?w=?(wTXTXw?yTXw?wTXTy+yTy)?w=XTXw+(XTX)Tw?(yTX)T?XTy=2XTXw?2XTy
其中XTXw=(XTX)TwXTXw=(XTX)Tw.令求導(dǎo)結(jié)果等于0，即可以求導(dǎo)問題的最小值。
2XTXw?2XTy=0w=(XTX)?1XTy
2XTXw?2XTy=0w=(XTX)?1XTy
??再利用最小二乘法的矩陣形式對前面的例子進行求解，用二階多項式擬合即兩個點(2,3),(5,8)(2,3),(5,8)。
表示輸入矩陣 XX和標(biāo)簽向量yy
X=[1125425]y=[38]T
X=[1241525]y=[38]T
計算XTXXTX
XTX=???272972913329133641???
XTX=[272972913329133641]
矩陣求逆，再做矩陣乘法運算
但 XTXXTX不可逆，故無唯一解。
??關(guān)于矩陣的逆是否存在，可以通過判斷矩陣的行列式是否為0（det(A)=?0det(A)=?0 來判斷，也可以通過初等行變換，觀察矩陣的行向量是否線性相關(guān)，在這個例子下，矩陣不可逆，故有無窮多解。但如果新增一個點(4,7)(4,7)，則就可以解了。
??其實這和數(shù)據(jù)集的點數(shù)和選擇的階數(shù)有關(guān)，如果點數(shù)小于階數(shù)則會出現(xiàn)無窮解的情況，如果點數(shù)等于階數(shù)，那么剛好有解可以完全擬合所有數(shù)據(jù)點，如果點數(shù)大于階數(shù)，則會求的近似解。
??那么對于點數(shù)小于階數(shù)的情況，如何求解？在python的多項式擬合函數(shù)中是可以擬合的，而且效果不錯，具體算法不是很了解，可以想辦法參考python的ployfit()函數(shù)的實現(xiàn)。
4. 擬合階數(shù)的選擇
?? 在前面的推導(dǎo)中，多項式的階數(shù)被固定了，那么實際場景下應(yīng)該如何選擇合適的階數(shù)MM呢？
一般會選擇階數(shù)MM小于點數(shù)NN
把訓(xùn)練數(shù)據(jù)分為訓(xùn)練集合驗證集，在訓(xùn)練集上，同時用不同的MM值訓(xùn)練多個模型，然后選擇在驗證集誤差最小的階數(shù)script type="math/tex" id="MathJax-Element-5573"M/script
求對 x 的偏導(dǎo)數(shù)，視 y 為常量，對 x 求導(dǎo)；
求對 y 的偏導(dǎo)數(shù)，視 x 為常量， 對 y 求導(dǎo)。
則：?f/?x = 4-2x, ?f/?y = -4-2y
偏導(dǎo)數(shù) f'x(x0,y0) 表示固定面上一點對 x 軸的切線斜率；偏導(dǎo)數(shù) f'y(x0,y0) 表示固定面上一點對 y 軸的切線斜率。
擴展資料：
將多元函數(shù)關(guān)于一個自變量求偏導(dǎo)數(shù)時，就將其余的自變量看成常數(shù)，此時求導(dǎo)方法與一元函數(shù)導(dǎo)數(shù)的求法是一樣的。
把 x 固定在 x0，讓 y 有增量 △y ，如果極限存在那么此極限稱為函數(shù) z=(x,y) 在 (x0,y0)處對 y 的偏導(dǎo)數(shù)。
偏函數(shù)是將所要承載的函數(shù)作為partial()函數(shù)的第一個參數(shù)，原函數(shù)的各個參數(shù)依次作為partial()函數(shù)后續(xù)的參數(shù)，除非使用關(guān)鍵字參數(shù)。
通過語言描述可能無法理解偏函數(shù)是怎么使用的，那么就舉一個常見的例子來說明。在這個例子里，我們實現(xiàn)了一個取余函數(shù)，對于整數(shù)100，取得對于不同數(shù)m的100%m的余數(shù)。
print（“字符串”），5/2和5//2的結(jié)果是不同的5/2為2.5,5//2為2.
python2需要導(dǎo)入from_future_import division執(zhí)行普通的除法。
1/2和1//2的結(jié)果0.5和0.
%號為取模運算。
乘方運算為2**3，-2**3和-（2**3）是等價的。
from sympy import*導(dǎo)入庫
x,y,z=symbols('x y z'),定義變量
init_printing(use_unicode=True)設(shè)置打印方式。
python的內(nèi)部常量有pi，
函數(shù)simplify，simplify(sin(x)**2 + cos(x)**2)化簡結(jié)果為1，
simplify((x**3 + x**2 - x - 1)/(x**2 + 2*x + 1))化簡結(jié)果為x-1?；嗁ゑR函數(shù)。simplify(gamma(x)/gamma(x - 2))得（x-2）(x-1)。
expand((x + 1)**2)展開多項式。
expand((x + 1)*(x - 2) - (x - 1)*x)
因式分解。factor(x**2*z + 4*x*y*z + 4*y**2*z)得到z*(x + 2*y)**2
from_future_import division
x,y,z,t=symbols('x y z t')定義變量，
k, m, n = symbols('k m n', integer=True)定義三個整數(shù)變量。
f, g, h = symbols('f g h', cls=Function)定義的類型為函數(shù)。
factor_list(x**2*z + 4*x*y*z + 4*y**2*z)得到一個列表，表示因式的冪，(1, [(z, 1), (x + 2*y, 2)])
expand((cos(x) + sin(x))**2)展開多項式。
expr = x*y + x - 3 + 2*x**2 - z*x**2 + x**3，collected_expr = collect(expr, x)將x合并。將x元素按階次整合。
collected_expr.coeff(x, 2)直接取出變量collected_expr的x的二次冪的系數(shù)。
cancel()is more efficient thanfactor().
cancel((x**2 + 2*x + 1)/(x**2 + x))
，expr = (x*y**2 - 2*x*y*z + x*z**2 + y**2 - 2*y*z + z**2)/(x**2 - 1)，cancel(expr)
expr = (4*x**3 + 21*x**2 + 10*x + 12)/(x**4 + 5*x**3 + 5*x**2 + 4*x)，apart(expr)
asin(1)
trigsimp(sin(x)**2 + cos(x)**2)三角函數(shù)表達式化簡，
trigsimp(sin(x)**4 - 2*cos(x)**2*sin(x)**2 + cos(x)**4)
trigsimp(sin(x)*tan(x)/sec(x))
trigsimp(cosh(x)**2 + sinh(x)**2)雙曲函數(shù)。
三角函數(shù)展開，expand_trig(sin(x + y))，acos(x)，cos(acos(x))，expand_trig(tan(2*x))
x, y = symbols('x y', positive=True)正數(shù)，a, b = symbols('a b', real=True)實數(shù)，z, t, c = symbols('z t c')定義變量的方法。
sqrt(x) == x**Rational(1, 2)判斷是否相等。
powsimp(x**a*x**b)冪函數(shù)的乘法，不同冪的乘法，必須先定義a和b。powsimp(x**a*y**a)相同冪的乘法。
powsimp(t**c*z**c)，注意，powsimp()refuses to do the simplification if it is not valid.
powsimp(t**c*z**c, force=True)這樣的話就可以得到化簡過的式子。聲明強制進行化簡。
(z*t)**2，sqrt(x*y)
第一個展開expand_power_exp(x**(a + b))，expand_power_base((x*y)**a)展開，
expand_power_base((z*t)**c, force=True)強制展開。
powdenest((x**a)**b)，powdenest((z**a)**b)，powdenest((z**a)**b, force=True)
ln(x)，x, y ,z= symbols('x y z', positive=True)，n = symbols('n', real=True)，
expand_log(log(x*y))展開為log(x) + log(y)，但是python3沒有。這是因為需要將x定義為positive。這是必須的，否則不會被展開。expand_log(log(x/y)),expand_log(log(x**n))
As withpowsimp()andpowdenest(),expand_log()has aforceoption that can be used to ignore assumptions。
expand_log(log(z**2), force=True),強制展開。
logcombine(log(x) + log(y))，logcombine(n*log(x))，logcombine(n*log(z), force=True)。
factorial(n)階乘，binomial(n, k)等于c（n，k），gamma(z)伽馬函數(shù)。
hyper([1, 2], [3], z)，
tan(x).rewrite(sin)得到用正弦表示的正切。factorial(x).rewrite(gamma)用伽馬函數(shù)重寫階乘。
expand_func(gamma(x + 3))得到，x*(x + 1)*(x + 2)*gamma(x)，
hyperexpand(hyper([1, 1], [2], z))，
combsimp(factorial(n)/factorial(n - 3))化簡，combsimp(binomial(n+1, k+1)/binomial(n, k))化簡。combsimp(gamma(x)*gamma(1 - x))
自定義函數(shù)
def list_to_frac(l):
expr = Integer(0)
for i in reversed(l[1:]):
expr += i
expr = 1/expr
return l[0] + expr
list_to_frac([x, y, z])結(jié)果為x + 1/z，這個結(jié)果是錯誤的。
syms = symbols('a0:5')，定義syms，得到的結(jié)果為(a0, a1, a2, a3, a4)。
這樣也可以a0, a1, a2, a3, a4 = syms， 可能是我的操作錯誤 。發(fā)現(xiàn)python和自動縮進有關(guān)，所以一定看好自動縮進的距離。list_to_frac([1, 2, 3, 4])結(jié)果為43/30。
使用cancel可以將生成的分式化簡，frac = cancel(frac)化簡為一個分?jǐn)?shù)線的分式。
(a0*a1*a2*a3*a4 + a0*a1*a2 + a0*a1*a4 + a0*a3*a4 + a0 + a2*a3*a4 + a2 + a4)/(a1*a2*a3*a4 + a1*a2 + a1*a4 + a3*a4 + 1)
a0, a1, a2, a3, a4 = syms定義a0到a4，frac = apart(frac, a0)可將a0提出來。frac=1/(frac-a0)將a0去掉取倒。frac = apart(frac, a1)提出a1。
help("modules"),模塊的含義，help("modules yourstr")模塊中包含的字符串的意思。，
help("topics")，import os.path + help("os.path")，help("list")，help("open")
# -*- coding: UTF-8 -*-聲明之后就可以在ide中使用中文注釋。
定義
l = list(symbols('a0:5'))定義列表得到[a0, a1, a2, a3, a4]
fromsympyimport*
x,y,z=symbols('x y z')
init_printing(use_unicode=True)
diff(cos(x),x)求導(dǎo)。diff(exp(x**2), x)，diff(x**4, x, x, x)和diff(x**4, x, 3)等價。
diff(expr, x, y, 2, z, 4)求出表達式的y的2階，z的4階，x的1階導(dǎo)數(shù)。和diff(expr, x, y, y, z, 4)等價。expr.diff(x, y, y, z, 4)一步到位。deriv = Derivative(expr, x, y, y, z, 4)求偏導(dǎo)。但是不顯示。之后用deriv.doit()即可顯示
integrate(cos(x), x)積分。定積分integrate(exp(-x), (x, 0, oo))無窮大用2個oo表示。integrate(exp(-x**2-y**2),(x,-oo,oo),(y,-oo,oo))二重積分。print(expr)print的使用。
expr = Integral(log(x)**2, x)，expr.doit()積分得到x*log(x)**2 - 2*x*log(x) + 2*x。
integ.doit()和integ = Integral((x**4 + x**2*exp(x) - x**2 - 2*x*exp(x) - 2*x -
exp(x))*exp(x)/((x - 1)**2*(x + 1)**2*(exp(x) + 1)), x)連用。
limit(sin(x)/x,x,0)，not-a-number表示nan算不出來，limit(expr, x, oo)，，expr = Limit((cos(x) - 1)/x, x, 0)，expr.doit()連用。左右極限limit(1/x, x, 0, '+')，limit(1/x, x, 0, '-')。。
Series Expansion級數(shù)展開。expr = exp(sin(x))，expr.series(x, 0, 4)得到1 + x + x**2/2 + O(x**4)，，x*O(1)得到O(x)，，expr.series(x, 0, 4).removeO()將無窮小移除。exp(x-6).series(x,x0=6)，，得到
-5 + (x - 6)**2/2 + (x - 6)**3/6 + (x - 6)**4/24 + (x - 6)**5/120 + x + O((x - 6)**6, (x, 6))最高到5階。
f=Function('f')定義函數(shù)變量和h=Symbol('h')和d2fdx2=f(x).diff(x,2)求2階，，as_finite_diff(dfdx)函數(shù)和as_finite_diff(d2fdx2,[-3*h,-h,2*h])，，x_list=[-3,1,2]和y_list=symbols('a b c')和apply_finite_diff(1,x_list,y_list,0)。
Eq(x, y)，，solveset(Eq(x**2, 1), x)解出來x，當(dāng)二式相等。和solveset(Eq(x**2 - 1, 0), x)等價。solveset(x**2 - 1, x)
solveset(x**2 - x, x)解，solveset(x - x, x, domain=S.Reals)解出來定義域。solveset(exp(x), x)? ? # No solution exists解出EmptySet()表示空集。
等式形式linsolve([x + y + z - 1, x + y + 2*z - 3 ], (x, y, z))和矩陣法linsolve(Matrix(([1, 1, 1, 1], [1, 1, 2, 3])), (x, y, z))得到{(-y - 1, y, 2)}
A*x = b 形式，M=Matrix(((1,1,1,1),(1,1,2,3)))，system=A,b=M[:,:-1],M[:,-1]，linsolve(system,x,y,z)，，solveset(x**3 - 6*x**2 + 9*x, x)解多項式。roots(x**3 - 6*x**2 + 9*x, x)，得出，{3: 2, 0: 1}，有2個3的重根，1個0根。solve([x*y - 1, x - 2], x, y)解出坐標(biāo)。
f, g = symbols('f g', cls=Function)函數(shù)的定義，解微分方程diffeq = Eq(f(x).diff(x, x) - 2*f(x).diff(x) + f(x), sin(x))再和dsolve(diffeq,f(x))結(jié)合。得到Eq(f(x), (C1 + C2*x)*exp(x) + cos(x)/2)，dsolve(f(x).diff(x)*(1 - sin(f(x))), f(x))解出來Eq(f(x) + cos(f(x)), C1)，，
Matrix([[1,-1],[3,4],[0,2]])，，Matrix([1, 2, 3])列表示。M=Matrix([[1,2,3],[3,2,1]])
N=Matrix([0,1,1])
M*N符合矩陣的乘法。M.shape顯示矩陣的行列數(shù)。
M.row(0)獲取M的第0行。M.col(-1)獲取倒數(shù)第一列。
M.col_del(0)刪掉第1列。M.row_del(1)刪除第二行，序列是從0開始的。M = M.row_insert(1, Matrix([[0, 4]]))插入第二行，，M = M.col_insert(0, Matrix([1, -2]))插入第一列。
M+N矩陣相加，M*N，3*M，M**2，M**-1，N**-1表示求逆。M.T求轉(zhuǎn)置。
eye(3)單位。zeros(2, 3)，0矩陣，ones(3, 2)全1，diag(1, 2, 3)對角矩陣。diag(-1, ones(2, 2), Matrix([5, 7, 5]))生成Matrix([
[-1, 0, 0, 0],
[ 0, 1, 1, 0],
[ 0, 1, 1, 0],
[ 0, 0, 0, 5],
[ 0, 0, 0, 7],
[ 0, 0, 0, 5]])矩陣。
Matrix([[1, 0, 1], [2, -1, 3], [4, 3, 2]])
一行一行顯示，，M.det()求行列式。M.rref()矩陣化簡。得到結(jié)果為Matrix([
[1, 0,? 1,? 3],
[0, 1, 2/3, 1/3],
[0, 0,? 0,? 0]]), [0, 1])。
M = Matrix([[1, 2, 3, 0, 0], [4, 10, 0, 0, 1]])，M.nullspace()
Columnspace
M.columnspace()和M = Matrix([[1, 2, 3, 0, 0], [4, 10, 0, 0, 1]])
M = Matrix([[3, -2,? 4, -2], [5,? 3, -3, -2], [5, -2,? 2, -2], [5, -2, -3,? 3]])和M.eigenvals()得到{3: 1, -2: 1, 5: 2}，，This means thatMhas eigenvalues -2, 3, and 5, and that the eigenvalues -2 and 3 have algebraic multiplicity 1 and that the eigenvalue 5 has algebraic multiplicity 2.
P, D = M.diagonalize()，P得Matrix([
[0, 1, 1,? 0],
[1, 1, 1, -1],
[1, 1, 1,? 0],
[1, 1, 0,? 1]])，，D為Matrix([
[-2, 0, 0, 0],
[ 0, 3, 0, 0],
[ 0, 0, 5, 0],
[ 0, 0, 0, 5]])
P*D*P**-1 == M返回為True。lamda = symbols('lamda')。
lamda = symbols('lamda')定義變量，p = M.charpoly(lamda)和factor(p)
expr = x**2 + x*y，srepr(expr)可以將表達式說明計算法則，"Add(Pow(Symbol('x'), Integer(2)), Mul(Symbol('x'), Symbol('y')))"。。
x = symbols('x')和x = Symbol('x')是一樣的。srepr(x**2)得到"Pow(Symbol('x'), Integer(2))"。Pow(x, 2)和Mul(x, y)得到x**2。x*y
type(2)得到class 'int'，type(sympify(2))得到class 'sympy.core.numbers.Integer'..srepr(x*y)得到"Mul(Symbol('x'), Symbol('y'))"。。。
Add(Pow(x, 2), Mul(x, y))得到"Add(Mul(Integer(-1), Pow(Symbol('x'), Integer(2))), Mul(Rational(1, 2), sin(Mul(Symbol('x'), Symbol('y')))), Pow(Symbol('y'), Integer(-1)))"。。Pow函數(shù)為冪次。
expr = Add(x, x)，expr.func。。Integer(2).func，class 'sympy.core.numbers.Integer',,Integer(0).func和Integer(-1).func，，，expr = 3*y**2*x和expr.func得到class 'sympy.core.mul.Mul',,expr.args將表達式分解為得到(3, x, y**2)，，expr.func(*expr.args)合并。expr == expr.func(*expr.args)返回True。expr.args[2]得到y(tǒng)**2，expr.args[1]得到x，expr.args[0]得到3.。
expr.args[2].args得到(y, 2)。。y.args得到空括號。Integer(2).args得到空括號。
from sympy import *
E**(I*pi)+1，可以看出，I和E，pi已將在sympy內(nèi)已定義。
x=Symbol('x')，，expand( E**(I*x) )不能展開，expand(exp(I*x),complex=True)可以展開，得到I*exp(-im(x))*sin(re(x)) + exp(-im(x))*cos(re(x))，，x=Symbol("x",real=True)將x定義為實數(shù)。再展開expand(exp(I*x),complex=True)得到。I*sin(x) + cos(x)。。
tmp = series(exp(I*x), x, 0, 10)和pprint(tmp)打印出來可讀性好，print(tmp)可讀性不好。。pprint將公式用更好看的格式打印出來，，pprint( series( cos(x), x, 0, 10) )
integrate(x*sin(x), x)，，定積分integrate(x*sin(x), (x, 0, 2*pi))。。
用雙重積分求解球的體積。
x, y, r = symbols('x,y,r')和2 * integrate(sqrt(r*r-x**2), (x, -r, r))計算球的體積。計算不來，是因為sympy不知道r是大于0的。r = symbols('r', positive=True)這樣定義r即可。circle_area=2*integrate(sqrt(r**2-x**2),(x,-r,r))得到。circle_area=circle_area.subs(r,sqrt(r**2-x**2))將r替換。
integrate(circle_area,(x,-r,r))再積分即可。
expression.sub([(x,y),(y,x)])又換到原來的狀況了。
expression.subs(x, y)，，將算式中的x替換成y。。
expression.subs({x:y,u:v}) : 使用字典進行多次替換。。
expression.subs([(x,y),(u,v)]) : 使用列表進行多次替換。。
最優(yōu)化
為什么要做最優(yōu)化呢？因為在生活中，人們總是希望幸福值或其它達到一個極值，比如做生意時希望成本最小，收入最大，所以在很多商業(yè)情境中，都會遇到求極值的情況。
函數(shù)求根
這里「函數(shù)的根」也稱「方程的根」，或「函數(shù)的零點」。
先把我們需要的包加載進來。import numpy as npimport scipy as spimport scipy.optimize as optimport matplotlib.pyplot as plt%matplotlib inline
函數(shù)求根和最優(yōu)化的關(guān)系？什么時候函數(shù)是最小值或最大值？
兩個問題一起回答：最優(yōu)化就是求函數(shù)的最小值或最大值，同時也是極值，在求一個函數(shù)最小值或最大值時，它所在的位置肯定是導(dǎo)數(shù)為 0 的位置，所以要求一個函數(shù)的極值，必然要先求導(dǎo)，使其為 0，所以函數(shù)求根就是為了得到最大值最小值。
scipy.optimize 有什么方法可以求根？
可以用 scipy.optimize 中的 bisect 或 brentq 求根。f = lambda x: np.cos(x) - x # 定義一個匿名函數(shù)x = np.linspace(-5, 5, 1000) # 先生成 1000 個 xy = f(x) # 對應(yīng)生成 1000 個 f(x)plt.plot(x, y); # 看一下這個函數(shù)長什么樣子plt.axhline(0, color='k'); # 畫一根橫線，位置在 y=0
opt.bisect(f, -5, 5) # 求取函數(shù)的根0.7390851332155535plt.plot(x, y)plt.axhline(0, color='k')plt.scatter([_], [0], c='r', s=100); # 這里的 [_] 表示上一個 Cell 中的結(jié)果，這里是 x 軸上的位置，0 是 y 上的位置
求根有兩種方法，除了上面介紹的 bisect，還有 brentq，后者比前者快很多。%timeit opt.bisect(f, -5, 5)%timeit opt.brentq(f, -5, 5)10000 loops, best of 3: 157 s per loopThe slowest run took 11.65 times longer than the fastest. This could mean that an intermediate result is being cached.10000 loops, best of 3: 35.9 s per loop
函數(shù)求最小化
求最小值就是一個最優(yōu)化問題。求最大值時只需對函數(shù)做一個轉(zhuǎn)換，比如加一個負(fù)號，或者取倒數(shù)，就可轉(zhuǎn)成求最小值問題。所以兩者是同一問題。
初始值對最優(yōu)化的影響是什么？
舉例來說，先定義個函數(shù)。f = lambda x: 1-np.sin(x)/xx = np.linspace(-20., 20., 1000)y = f(x)
當(dāng)初始值為 3 值，使用 minimize 函數(shù)找到最小值。minimize 函數(shù)是在新版的 scipy 里，取代了以前的很多最優(yōu)化函數(shù)，是個通用的接口，背后是很多方法在支撐。x0 = 3xmin = opt.minimize(f, x0).x # x0 是起始點，起始點最好離真正的最小值點不要太遠plt.plot(x, y)plt.scatter(x0, f(x0), marker='o', s=300); # 起始點畫出來，用圓圈表示plt.scatter(xmin, f(xmin), marker='v', s=300); # 最小值點畫出來，用三角表示plt.xlim(-20, 20);
初始值為 3 時，成功找到最小值。
現(xiàn)在來看看初始值為 10 時，找到的最小值點。x0 = 10xmin = opt.minimize(f, x0).xplt.plot(x, y)plt.scatter(x0, f(x0), marker='o', s=300)plt.scatter(xmin, f(xmin), marker='v', s=300)plt.xlim(-20, 20);
由上圖可見，當(dāng)初始值為 10 時，函數(shù)找到的是局部最小值點，可見 minimize 的默認(rèn)算法對起始點的依賴性。
那么怎么才能不管初始值在哪個位置，都能找到全局最小值點呢？
如何找到全局最優(yōu)點？
可以使用 basinhopping 函數(shù)找到全局最優(yōu)點，相關(guān)背后算法，可以看幫助文件，有提供論文的索引和出處。
我們設(shè)初始值為 10 看是否能找到全局最小值點。x0 = 10from scipy.optimize import basinhoppingxmin = basinhopping(f,x0,stepsize = 5).xplt.plot(x, y);plt.scatter(x0, f(x0), marker='o', s=300);plt.scatter(xmin, f(xmin), marker='v', s=300);plt.xlim(-20, 20);
當(dāng)起始點在比較遠的位置，依然成功找到了全局最小值點。
如何求多元函數(shù)最小值？
以二元函數(shù)為例，使用 minimize 求對應(yīng)的最小值。def g(X): x,y = X return (x-1)**4 + 5 * (y-1)**2 - 2*x*yX_opt = opt.minimize(g, (8, 3)).x # (8,3) 是起始點print X_opt[ 1.88292611 1.37658521]fig, ax = plt.subplots(figsize=(6, 4)) # 定義畫布和圖形x_ = y_ = np.linspace(-1, 4, 100)X, Y = np.meshgrid(x_, y_)c = ax.contour(X, Y, g((X, Y)), 50) # 等高線圖ax.plot(X_opt[0], X_opt[1], 'r*', markersize=15) # 最小點的位置是個元組ax.set_xlabel(r"$x_1$", fontsize=18)ax.set_ylabel(r"$x_2$", fontsize=18)plt.colorbar(c, ax=ax) # colorbar 表示顏色越深，高度越高fig.tight_layout()
畫3D 圖。from mpl_toolkits.mplot3d import Axes3Dfrom matplotlib import cmfig = plt.figure()ax = fig.gca(projection='3d')x_ = y_ = np.linspace(-1, 4, 100)X, Y = np.meshgrid(x_, y_)surf = ax.plot_surface(X, Y, g((X,Y)), rstride=1, cstride=1, cmap=cm.coolwarm, linewidth=0, antialiased=False)cset = ax.contour(X, Y, g((X,Y)), zdir='z',offset=-5, cmap=cm.coolwarm)fig.colorbar(surf, shrink=0.5, aspect=5);
曲線擬合
曲線擬合和最優(yōu)化有什么關(guān)系？
曲線擬合的問題是，給定一組數(shù)據(jù)，它可能是沿著一條線散布的，這時要找到一條最優(yōu)的曲線來擬合這些數(shù)據(jù)，也就是要找到最好的線來代表這些點，這里的最優(yōu)是指這些點和線之間的距離是最小的，這就是為什么要用最優(yōu)化問題來解決曲線擬合問題。
舉例說明，給一些點，找到一條線，來擬合這些點。
先給定一些點：N = 50 # 點的個數(shù)m_true = 2 # 斜率b_true = -1 # 截距dy = 2.0 # 誤差np.random.seed(0)xdata = 10 * np.random.random(N) # 50 個 x，服從均勻分布ydata = np.random.normal(b_true + m_true * xdata, dy) # dy 是標(biāo)準(zhǔn)差plt.errorbar(xdata, ydata, dy, fmt='.k', ecolor='lightgray');
上面的點整體上呈現(xiàn)一個線性關(guān)系，要找到一條斜線來代表這些點，這就是經(jīng)典的一元線性回歸。目標(biāo)就是找到最好的線，使點和線的距離最短。要優(yōu)化的函數(shù)是點和線之間的距離，使其最小。點是確定的，而線是可變的，線是由參數(shù)值，斜率和截距決定的，這里就是要通過優(yōu)化距離找到最優(yōu)的斜率和截距。
點和線的距離定義如下：def chi2(theta, x, y): return np.sum(((y - theta[0] - theta[1] * x)) ** 2)
上式就是誤差平方和。
誤差平方和是什么？有什么作用？
誤差平方和公式為：
誤差平方和大，表示真實的點和預(yù)測的線之間距離太遠，說明擬合得不好，最好的線，應(yīng)該是使誤差平方和最小，即最優(yōu)的擬合線，這里是條直線。
誤差平方和就是要最小化的目標(biāo)函數(shù)。
找到最優(yōu)的函數(shù)，即斜率和截距。theta_guess = [0, 1] # 初始值theta_best = opt.minimize(chi2, theta_guess, args=(xdata, ydata)).xprint(theta_best)[-1.01442005 1.93854656]
上面兩個輸出即是預(yù)測的直線斜率和截距，我們是根據(jù)點來反推直線的斜率和截距，那么真實的斜率和截距是多少呢？-1 和 2，很接近了，差的一點是因為有噪音的引入。xfit = np.linspace(0, 10)yfit = theta_best[0] + theta_best[1] * xfitplt.errorbar(xdata, ydata, dy, fmt='.k', ecolor='lightgray');plt.plot(xfit, yfit, '-k');
最小二乘（Least Square）是什么？
上面用的是 minimize 方法，這個問題的目標(biāo)函數(shù)是誤差平方和，這就又有一個特定的解法，即最小二乘。
最小二乘的思想就是要使得觀測點和估計點的距離的平方和達到最小，這里的“二乘”指的是用平方來度量觀測點與估計點的遠近（在古漢語中“平方”稱為“二乘”），“最小”指的是參數(shù)的估計值要保證各個觀測點與估計點的距離的平方和達到最小。
關(guān)于最小二乘估計的計算，涉及更多的數(shù)學(xué)知識，這里不想詳述，其一般的過程是用目標(biāo)函數(shù)對各參數(shù)求偏導(dǎo)數(shù)，并令其等于 0，得到一個線性方程組。具體推導(dǎo)過程可參考斯坦福機器學(xué)習(xí)講義 第 7 頁。def deviations(theta, x, y): return (y - theta[0] - theta[1] * x)theta_best, ier = opt.leastsq(deviations, theta_guess, args=(xdata, ydata))print(theta_best)[-1.01442016 1.93854659]
最小二乘 leastsq 的結(jié)果跟 minimize 結(jié)果一樣。注意 leastsq 的第一個參數(shù)不再是誤差平方和 chi2，而是誤差本身 deviations，即沒有平方，也沒有和。yfit = theta_best[0] + theta_best[1] * xfitplt.errorbar(xdata, ydata, dy, fmt='.k', ecolor='lightgray');plt.plot(xfit, yfit, '-k');
非線性最小二乘
上面是給一些點，擬合一條直線，擬合一條曲線也是一樣的。def f(x, beta0, beta1, beta2): # 首先定義一個非線性函數(shù)，有 3 個參數(shù) return beta0 + beta1 * np.exp(-beta2 * x**2)beta = (0.25, 0.75, 0.5) # 先猜 3 個 betaxdata = np.linspace(0, 5, 50)y = f(xdata, *beta)ydata = y + 0.05 * np.random.randn(len(xdata)) # 給 y 加噪音def g(beta): return ydata - f(xdata, *beta) # 真實 y 和 預(yù)測值的差，求最優(yōu)曲線時要用到beta_start = (1, 1, 1)beta_opt, beta_cov = opt.leastsq(g, beta_start)print beta_opt # 求到的 3 個最優(yōu)的 beta 值[ 0.25525709 0.74270226 0.54966466]
拿估計的 beta_opt 值跟真實的 beta = (0.25, 0.75, 0.5) 值比較，差不多。fig, ax = plt.subplots()ax.scatter(xdata, ydata) # 畫點ax.plot(xdata, y, 'r', lw=2) # 真實值的線ax.plot(xdata, f(xdata, *beta_opt), 'b', lw=2) # 擬合的線ax.set_xlim(0, 5)ax.set_xlabel(r"$x$", fontsize=18)ax.set_ylabel(r"$f(x, \beta)$", fontsize=18)fig.tight_layout()
除了使用最小二乘，還可以使用曲線擬合的方法，得到的結(jié)果是一樣的。beta_opt, beta_cov = opt.curve_fit(f, xdata, ydata)print beta_opt[ 0.25525709 0.74270226 0.54966466]
有約束的最小化
有約束的最小化是指，要求函數(shù)最小化之外，還要滿足約束條件，舉例說明。
邊界約束def f(X): x, y = X return (x-1)**2 + (y-1)**2 # 這是一個碗狀的函數(shù)x_opt = opt.minimize(f, (0, 0), method='BFGS').x # 無約束最優(yōu)化
假設(shè)有約束條件，x 和 y 要在一定的范圍內(nèi)，如 x 在 2 到 3 之間，y 在 0 和 2 之間。bnd_x1, bnd_x2 = (2, 3), (0, 2) # 對自變量的約束x_cons_opt = opt.minimize(f, np.array([0, 0]), method='L-BFGS-B', bounds=[bnd_x1, bnd_x2]).x # bounds 矩形約束fig, ax = plt.subplots(figsize=(6, 4))x_ = y_ = np.linspace(-1, 3, 100)X, Y = np.meshgrid(x_, y_)c = ax.contour(X, Y, f((X,Y)), 50)ax.plot(x_opt[0], x_opt[1], 'b*', markersize=15) # 沒有約束下的最小值，藍色五角星ax.plot(x_cons_opt[0], x_cons_opt[1], 'r*', markersize=15) # 有約束下的最小值，紅色星星bound_rect = plt.Rectangle((bnd_x1[0], bnd_x2[0]), bnd_x1[1] - bnd_x1[0], bnd_x2[1] - bnd_x2[0], facecolor="grey")ax.add_patch(bound_rect)ax.set_xlabel(r"$x_1$", fontsize=18)ax.set_ylabel(r"$x_2$", fontsize=18)plt.colorbar(c, ax=ax)fig.tight_layout()
不等式約束
介紹下相關(guān)理論，先來看下存在等式約束的極值問題求法，比如下面的優(yōu)化問題。
目標(biāo)函數(shù)是 f(w)，下面是等式約束，通常解法是引入拉格朗日算子，這里使用 ββ 來表示算子，得到拉格朗日公式為
l 是等式約束的個數(shù)。
然后分別對 w 和ββ 求偏導(dǎo)，使得偏導(dǎo)數(shù)等于 0，然后解出 w 和βiβi，至于為什么引入拉格朗日算子可以求出極值，原因是 f(w) 的 dw 變化方向受其他不等式的約束，dw的變化方向與f(w)的梯度垂直時才能獲得極值，而且在極值處，f(w) 的梯度與其他等式梯度的線性組合平行，因此他們之間存在線性關(guān)系。（參考《最優(yōu)化與KKT條件》）
對于不等式約束的極值問題
常常利用拉格朗日對偶性將原始問題轉(zhuǎn)換為對偶問題，通過解對偶問題而得到原始問題的解。該方法應(yīng)用在許多統(tǒng)計學(xué)習(xí)方法中。有興趣的可以參閱相關(guān)資料，這里不再贅述。def f(X): return (X[0] - 1)**2 + (X[1] - 1)**2def g(X): return X[1] - 1.75 - (X[0] - 0.75)**4x_opt = opt.minimize(f, (0, 0), method='BFGS').xconstraints = [dict(type='ineq', fun=g)] # 約束采用字典定義，約束方式為不等式約束，邊界用 g 表示x_cons_opt = opt.minimize(f, (0, 0), method='SLSQP', constraints=constraints).xfig, ax = plt.subplots(figsize=(6, 4))x_ = y_ = np.linspace(-1, 3, 100)X, Y = np.meshgrid(x_, y_)c = ax.contour(X, Y, f((X, Y)), 50)ax.plot(x_opt[0], x_opt[1], 'b*', markersize=15) # 藍色星星，沒有約束下的最小值ax.plot(x_, 1.75 + (x_-0.75)**4, '', markersize=15)ax.fill_between(x_, 1.75 + (x_-0.75)**4, 3, color="grey")ax.plot(x_cons_opt[0], x_cons_opt[1], 'r*', markersize=15) # 在區(qū)域約束下的最小值ax.set_ylim(-1, 3)ax.set_xlabel(r"$x_0$", fontsize=18)ax.set_ylabel(r"$x_1$", fontsize=18)plt.colorbar(c, ax=ax)fig.tight_layout()
scipy.optimize.minimize 中包括了多種最優(yōu)化算法，每種算法使用范圍不同，詳細(xì)參考官方文檔。
重慶分公司
重慶分公司
