一、SciPy概述
前篇已經(jīng)大致介紹了NumPy，接下來(lái)讓我們看看SciPy能做些什么。NumPy替我們搞定了向量和矩陣的相關(guān)操作，基本上算是一個(gè)高級(jí)的科學(xué)計(jì)算器。SciPy基于NumPy提供了更為豐富和高級(jí)的功能擴(kuò)展，在統(tǒng)計(jì)、優(yōu)化、插值、數(shù)值積分、時(shí)頻轉(zhuǎn)換等方面提供了大量的可用函數(shù)，基本覆蓋了基礎(chǔ)科學(xué)計(jì)算相關(guān)的問(wèn)題。

在量化分析中，運(yùn)用最廣泛的是統(tǒng)計(jì)和優(yōu)化的相關(guān)技術(shù)，本篇重點(diǎn)介紹SciPy中的統(tǒng)計(jì)和優(yōu)化模塊，其他模塊在隨后系列文章中用到時(shí)再做詳述。
本篇會(huì)涉及到一些矩陣代數(shù)，如若感覺(jué)不適，可考慮跳過(guò)第三部分或者在理解時(shí)簡(jiǎn)單采用一維的標(biāo)量代替高維的向量。
首先還是導(dǎo)入相關(guān)的模塊，我們使用的是SciPy里面的統(tǒng)計(jì)和優(yōu)化部分：
In [1]:

import numpy as np
import scipy.stats as stats
import scipy.optimize as opt
二、統(tǒng)計(jì)部分?
2.1 生成隨機(jī)數(shù)?
我們從生成隨機(jī)數(shù)開(kāi)始，這樣方便后面的介紹。生成n個(gè)隨機(jī)數(shù)可用rv_continuous.rvs(size=n)或rv_discrete.rvs(size=n)，其中rv_continuous表示連續(xù)型的隨機(jī)分布，如均勻分布（uniform）、正態(tài)分布（norm）、貝塔分布（beta）等；rv_discrete表示離散型的隨機(jī)分布，如伯努利分布（bernoulli）、幾何分布（geom）、泊松分布（poisson）等。我們生成10個(gè)[0, 1]區(qū)間上的隨機(jī)數(shù)和10個(gè)服從參數(shù)$a = 4$，$b = 2$的貝塔分布隨機(jī)數(shù)：
In [2]:

rv_unif = stats.uniform.rvs(size=10)
print rv_unif
rv_beta = stats.beta.rvs(size=10, a=4, b=2)
print rv_beta

 

[ 0.20630272  0.25929204  0.16859206  0.92573462  0.16383319  0.3475617
  0.83792048  0.79574153  0.37945051  0.23439682]
[ 0.71216492  0.85688464  0.70310131  0.3783662   0.69507561  0.78626586
  0.54529967  0.4261079   0.26646767  0.8519046 ]
在每個(gè)隨機(jī)分布的生成函數(shù)里，都內(nèi)置了默認(rèn)的參數(shù)，如均勻分布的上下界默認(rèn)是0和1。可是一旦需要修改這些參數(shù)，每次生成隨機(jī)都要敲這么老長(zhǎng)一串有點(diǎn)麻煩，能不能簡(jiǎn)單點(diǎn)？SciPy里頭有一個(gè)Freezing的功能，可以提供簡(jiǎn)便版本的命令。SciPy.stats支持定義出某個(gè)具體的分布的對(duì)象，我們可以做如下的定義，讓beta直接指代具體參數(shù)$a = 4$和$b = 2$的貝塔分布。為讓結(jié)果具有可比性，這里指定了隨機(jī)數(shù)的生成種子，由NumPy提供。
In [3]:

np.random.seed(seed=2015)
rv_beta = stats.beta.rvs(size=10, a=4, b=2)
print "method 1:"
print rv_beta

np.random.seed(seed=2015)
beta = stats.beta(a=4, b=2)
print "method 2:"
print beta.rvs(size=10)

method 1:
[ 0.43857338  0.9411551   0.75116671  0.92002864  0.62030521  0.56585548
  0.41843548  0.5953096   0.88983036  0.94675351]
method 2:
[ 0.43857338  0.9411551   0.75116671  0.92002864  0.62030521  0.56585548
  0.41843548  0.5953096   0.88983036  0.94675351]
2.2 假設(shè)檢驗(yàn)?
好了，現(xiàn)在我們生成一組數(shù)據(jù)，并查看相關(guān)的統(tǒng)計(jì)量（相關(guān)分布的參數(shù)可以在這里查到：http://docs.scipy.org/doc/scipy/reference/stats.html）：
In [4]:

norm_dist = stats.norm(loc=0.5, scale=2)
n = 200
dat = norm_dist.rvs(size=n)
print "mean of data is: " + str(np.mean(dat))
print "median of data is: " + str(np.median(dat))
print "standard deviation of data is: " + str(np.std(dat))

mean of data is: 0.705195138069
median of data is: 0.658167882933
standard deviation of data is: 2.08967006905
假設(shè)這個(gè)數(shù)據(jù)是我們獲取到的實(shí)際的某些數(shù)據(jù)，如股票日漲跌幅，我們對(duì)數(shù)據(jù)進(jìn)行簡(jiǎn)單的分析。最簡(jiǎn)單的是檢驗(yàn)這一組數(shù)據(jù)是否服從假設(shè)的分布，如正態(tài)分布。這個(gè)問(wèn)題是典型的單樣本假設(shè)檢驗(yàn)問(wèn)題，最為常見(jiàn)的解決方案是采用K-S檢驗(yàn)（ Kolmogorov-Smirnov test）。單樣本K-S檢驗(yàn)的原假設(shè)是給定的數(shù)據(jù)來(lái)自和原假設(shè)分布相同的分布，在SciPy中提供了kstest函數(shù)，參數(shù)分別是數(shù)據(jù)、擬檢驗(yàn)的分布名稱和對(duì)應(yīng)的參數(shù)：
In [5]:

mu = np.mean(dat)
sigma = np.std(dat)
stat_val, p_val = stats.kstest(dat, 'norm', (mu, sigma))
print 'KS-statistic D = %6.3f p-value = %6.4f' % (stat_val, p_val)

KS-statistic D =  0.045 p-value = 0.8195

假設(shè)檢驗(yàn)的$p$-value值很大（在原假設(shè)下，$p$-value是服從[0, 1]區(qū)間上的均勻分布的隨機(jī)變量，可參考http://en.wikipedia.org/wiki/P-value ），因此我們接受原假設(shè)，即該數(shù)據(jù)通過(guò)了正態(tài)性的檢驗(yàn)。在正態(tài)性的前提下，我們可進(jìn)一步檢驗(yàn)這組數(shù)據(jù)的均值是不是0。典型的方法是$t$檢驗(yàn)（$t$-test），其中單樣本的$t$檢驗(yàn)函數(shù)為ttest_1samp：
In [6]:

stat_val, p_val = stats.ttest_1samp(dat, 0)
print 'One-sample t-statistic D = %6.3f, p-value = %6.4f' % (stat_val, p_val)

One-sample t-statistic D =  4.761, p-value = 0.0000

我們看到$p$-value$ < 0.05$，即給定顯著性水平0.05的前提下，我們應(yīng)拒絕原假設(shè)：數(shù)據(jù)的均值為0。我們?cè)偕梢唤M數(shù)據(jù)，嘗試一下雙樣本的$t$檢驗(yàn)（ttest_ind）：
In [7]:

norm_dist2 = stats.norm(loc=-0.2, scale=1.2)
dat2 = norm_dist2.rvs(size=n/2)
stat_val, p_val = stats.ttest_ind(dat, dat2, equal_var=False)
print 'Two-sample t-statistic D = %6.3f, p-value = %6.4f' % (stat_val, p_val)

Two-sample t-statistic D =  5.565, p-value = 0.0000
注意，這里我們生成的第二組數(shù)據(jù)樣本大小、方差和第一組均不相等，在運(yùn)用$t$檢驗(yàn)時(shí)需要使用Welch's $t$-test，即指定ttest_ind中的equal_var=False。我們同樣得到了比較小的$p$-value$，在顯著性水平0.05的前提下拒絕原假設(shè)，即認(rèn)為兩組數(shù)據(jù)均值不等。

stats還提供其他大量的假設(shè)檢驗(yàn)函數(shù)，如bartlett和levene用于檢驗(yàn)方差是否相等；anderson_ksamp用于進(jìn)行Anderson-Darling的K-樣本檢驗(yàn)等。

2.3 其他函數(shù)?

有時(shí)需要知道某數(shù)值在一個(gè)分布中的分位，或者給定了一個(gè)分布，求某分位上的數(shù)值。這可以通過(guò)cdf和ppf函數(shù)完成：
In [8]:

g_dist = stats.gamma(a=2)
print "quantiles of 2, 4 and 5:"
print g_dist.cdf([2, 4, 5])
print "Values of 25%, 50% and 90%:"

print g_dist.pdf([0.25, 0.5, 0.95])

quantiles of 2, 4 and 5:
[ 0.59399415  0.90842181  0.95957232]
Values of 25%, 50% and 90%:
[ 0.1947002   0.30326533  0.36740397]

對(duì)于一個(gè)給定的分布，可以用moment很方便的查看分布的矩信息，例如我們查看$N(0, 1)$的六階原點(diǎn)矩：
In [9]:

stats.norm.moment(6, loc=0, scale=1)

Out[9]:

15.000000000895332

describe函數(shù)提供對(duì)數(shù)據(jù)集的統(tǒng)計(jì)描述分析，包括數(shù)據(jù)樣本大小，極值，均值，方差，偏度和峰度：
In [10]:

norm_dist = stats.norm(loc=0, scale=1.8)
dat = norm_dist.rvs(size=100)
info = stats.describe(dat)
print "Data size is: " + str(info[0])
print "Minimum value is: " + str(info[1][0])
print "Maximum value is: " + str(info[1][1])
print "Arithmetic mean is: " + str(info[2])
print "Unbiased variance is: " + str(info[3])
print "Biased skewness is: " + str(info[4])

print "Biased kurtosis is: " + str(info[5])

Data size is: 100
Minimum value is: -4.12414564687
Maximum value is: 4.82577602489
Arithmetic mean is: 0.0962913592209
Unbiased variance is: 2.88719292463
Biased skewness is: -0.00256548794681
Biased kurtosis is: -0.317463421177

當(dāng)我們知道一組數(shù)據(jù)服從某些分布的時(shí)候，可以調(diào)用fit函數(shù)來(lái)得到對(duì)應(yīng)分布參數(shù)的極大似然估計(jì)（MLE, maximum-likelihood estimation）。以下代碼示例了假設(shè)數(shù)據(jù)服從正態(tài)分布，用極大似然估計(jì)分布參數(shù)：
In [11]:

norm_dist = stats.norm(loc=0, scale=1.8)
dat = norm_dist.rvs(size=100)
mu, sigma = stats.norm.fit(dat)
print "MLE of data mean:" + str(mu)
print "MLE of data standard deviation:" + str(sigma)

MLE of data mean:-0.249880829912
MLE of data standard deviation:1.89195303507

pearsonr和spearmanr可以計(jì)算Pearson和Spearman相關(guān)系數(shù)，這兩個(gè)相關(guān)系數(shù)度量了兩組數(shù)據(jù)的相互線性關(guān)聯(lián)程度：
In [12]:

norm_dist = stats.norm()
dat1 = norm_dist.rvs(size=100)
exp_dist = stats.expon()
dat2 = exp_dist.rvs(size=100)
cor, pval = stats.pearsonr(dat1, dat2)
print "Pearson correlation coefficient: " + str(cor)
cor, pval = stats.pearsonr(dat1, dat2)
print "Spearman's rank correlation coefficient: " + str(cor)

Pearson correlation coefficient: -0.0262911931014
Spearman's rank correlation coefficient: -0.0262911931014

其中的$p$-value表示原假設(shè)（兩組數(shù)據(jù)不相關(guān)）下，相關(guān)系數(shù)的顯著性。

最后，在分析金融數(shù)據(jù)中使用頻繁的線性回歸在SciPy中也有提供，我們來(lái)看一個(gè)例子：
In [13]:

x = stats.chi2.rvs(3, size=50)
y = 2.5 + 1.2 * x + stats.norm.rvs(size=50, loc=0, scale=1.5)
slope, intercept, r_value, p_value, std_err = stats.linregress(x, y)
print "Slope of fitted model is:" , slope
print "Intercept of fitted model is:", intercept
print "R-squared:", r_value**2

Slope of fitted model is: 1.44515601191
Intercept of fitted model is: 1.91080684516
R-squared: 0.798786910173
在前面的鏈接中，可以查到大部分stat中的函數(shù)，本節(jié)權(quán)作簡(jiǎn)單介紹，挖掘更多功能的最好方法還是直接讀原始的文檔。另外，StatsModels（http://statsmodels.sourceforge.net ）模塊提供了更為專業(yè)，更多的統(tǒng)計(jì)相關(guān)函數(shù)。若在SciPy沒(méi)有滿足需求，可以采用StatsModels。

三、優(yōu)化部分?

優(yōu)化問(wèn)題在投資中可謂是根本問(wèn)題，如果手上有眾多可選的策略，應(yīng)如何從中選擇一個(gè)“最好”的策略進(jìn)行投資呢？這時(shí)就需要用到一些優(yōu)化技術(shù)針對(duì)給定的指標(biāo)進(jìn)行尋優(yōu)。隨著越來(lái)越多金融數(shù)據(jù)的出現(xiàn)，機(jī)器學(xué)習(xí)逐漸應(yīng)用在投資領(lǐng)域，在機(jī)器學(xué)習(xí)中，優(yōu)化也是十分重要的一個(gè)部分。以下介紹一些常見(jiàn)的優(yōu)化方法，雖然例子是人工生成的，不直接應(yīng)用于實(shí)際金融數(shù)據(jù)，我們希望讀者在后面遇到優(yōu)化問(wèn)題時(shí)，能夠從這些簡(jiǎn)單例子迅速上手解決。
 
3.1 無(wú)約束優(yōu)化問(wèn)題?

所謂的無(wú)約束優(yōu)化問(wèn)題指的是一個(gè)優(yōu)化問(wèn)題的尋優(yōu)可行集合是目標(biāo)函數(shù)自變量的定義域，即沒(méi)有外部的限制條件。例如，求解優(yōu)化問(wèn)題 [
minimizef(x)=x2?4.8x+1.2
] 就是一個(gè)無(wú)約束優(yōu)化問(wèn)題，而求解 [
minimizef(x)=x2?4.8x+1.2 subject tox≥0

]

則是一個(gè)帶約束的優(yōu)化問(wèn)題。更進(jìn)一步，我們假設(shè)考慮的問(wèn)題全部是凸優(yōu)化問(wèn)題，即目標(biāo)函數(shù)是凸函數(shù)，其自變量的可行集是凸集。（詳細(xì)定義可參考斯坦福大學(xué)Stephen Boyd教授的教材convex optimization，下載鏈接：http://stanford.edu/~boyd/cvxbook ）

我們以Rosenbrock函數(shù) [ f(\mathbf{x}) = \sum{i=1}^{N-1} 100 (x_i - x{i-1}^2)^2 + (1 - x_{i-1})^2 ] 作為尋優(yōu)的目標(biāo)函數(shù)來(lái)簡(jiǎn)要介紹在SciPy中使用優(yōu)化模塊scipy.optimize。

首先需要定義一下這個(gè)Rosenbrock函數(shù)：
In [14]:

def rosen(x):
    """The Rosenbrock function"""
    return sum(100.0*(x[1:]-x[:-1]**2.0)**2.0 + (1-x[:-1])**2.0)

3.1.1 Nelder-Mead單純形法?

單純形法是運(yùn)籌學(xué)中介紹的求解線性規(guī)劃問(wèn)題的通用方法，這里的Nelder-Mead單純形法與其并不相同，只是用到單純形的概念。設(shè)定起始點(diǎn)$\mathbf{x}_0 = (1.3, 0.7, 0.8, 1.9, 1.2)$，并進(jìn)行最小化的尋優(yōu)。這里‘xtol’表示迭代收斂的容忍誤差上界：
In [15]:

x_0 = np.array([0.5, 1.6, 1.1, 0.8, 1.2])
res = opt.minimize(rosen, x_0, method='nelder-mead', options={'xtol': 1e-8, 'disp': True})
print "Result of minimizing Rosenbrock function via Nelder-Mead Simplex algorithm:"
print res

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 436
         Function evaluations: 706
Result of minimizing Rosenbrock function via Nelder-Mead Simplex algorithm:
  status: 0
    nfev: 706
 success: True
     fun: 1.6614969876635003e-17
       x: array([ 1.,  1.,  1.,  1.,  1.])
 message: 'Optimization terminated successfully.'
     nit: 436

Rosenbrock函數(shù)的性質(zhì)比較好，簡(jiǎn)單的優(yōu)化方法就可以處理了，還可以在minimize中使用method='powell'來(lái)指定使用Powell's method。這兩種簡(jiǎn)單的方法并不使用函數(shù)的梯度，在略微復(fù)雜的情形下收斂速度比較慢，下面讓我們來(lái)看一下用到函數(shù)梯度進(jìn)行尋優(yōu)的方法。
 
3.1.2 Broyden-Fletcher-Goldfarb-Shanno法?

Broyden-Fletcher-Goldfarb-Shanno（BFGS）法用到了梯度信息，首先求一下Rosenbrock函數(shù)的梯度：

[ \begin{split} \frac{\partial f}{\partial xj} &= \sum{i=1}^N 200(xi - x{i-1}^2)(\delta{i,j} - 2x{i-1}\delta{i-1,j}) -2(1 - x{i-1})\delta_{i-1,j} \ &= 200(xj - x{j-1}^2) - 400xj(x{j+1} - x_j^2) - 2(1 - x_j) \end{split}] 其中當(dāng)$i=j$時(shí)，$\delta_{i,j} = 1$，否則$\delta_{i,j} = 0$。

邊界的梯度是特例，有如下形式： [ \begin{split} \frac{\partial f}{\partial x_0} &= -400x_0(x_1 - x_0^2) - 2(1 - x_0), \ \frac{\partial f}{\partial x{N-1}} &= 200(x{N-1} - x_{N-2}^2) \end{split}]

我們可以如下定義梯度向量的計(jì)算函數(shù)了：
In [16]:

def rosen_der(x):
    xm = x[1:-1]
    xm_m1 = x[:-2]
    xm_p1 = x[2:]
    der = np.zeros_like(x)
    der[1:-1] = 200*(xm-xm_m1**2) - 400*(xm_p1 - xm**2)*xm - 2*(1-xm)
    der[0] = -400*x[0]*(x[1]-x[0]**2) - 2*(1-x[0])
    der[-1] = 200*(x[-1]-x[-2]**2)
    return der
梯度信息的引入在minimize函數(shù)中通過(guò)參數(shù)jac指定：
In [17]:

res = opt.minimize(rosen, x_0, method='BFGS', jac=rosen_der, options={'disp': True})
print "Result of minimizing Rosenbrock function via Broyden-Fletcher-Goldfarb-Shanno algorithm:"
print res

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 52
         Function evaluations: 63
         Gradient evaluations: 63
Result of minimizing Rosenbrock function via Broyden-Fletcher-Goldfarb-Shanno algorithm:
   status: 0
  success: True
     njev: 63
     nfev: 63
 hess_inv: array([[ 0.00726515,  0.01195827,  0.0225785 ,  0.04460906,  0.08923649],
       [ 0.01195827,  0.02417936,  0.04591135,  0.09086889,  0.18165604],
       [ 0.0225785 ,  0.04591135,  0.09208689,  0.18237695,  0.36445491],
       [ 0.04460906,  0.09086889,  0.18237695,  0.36609277,  0.73152922],
       [ 0.08923649,  0.18165604,  0.36445491,  0.73152922,  1.46680958]])
      fun: 3.179561068096293e-14
        x: array([ 1.        ,  0.99999998,  0.99999996,  0.99999992,  0.99999983])
  message: 'Optimization terminated successfully.'
      jac: array([  4.47207141e-06,   1.30357917e-06,  -1.86454207e-07,
        -2.00564982e-06,   4.98799446e-07])

3.1.3 牛頓共軛梯度法（Newton-Conjugate-Gradient algorithm）?

用到梯度的方法還有牛頓法，牛頓法是收斂速度最快的方法，其缺點(diǎn)在于要求Hessian矩陣（二階導(dǎo)數(shù)矩陣）。牛頓法大致的思路是采用泰勒展開(kāi)的二階近似： [ f(\mathbf{x}) \approx f(\mathbf{x}_0) + \nabla f(\mathbf{x}_0)(\mathbf{x} - \mathbf{x}_0) + \frac{1}{2}(\mathbf{x} - \mathbf{x}_0)^T\mathbf{H}(\mathbf{x}_0)(\mathbf{x} - \mathbf{x}_0) ] 其中$\mathbf{H}(\mathbf{x}_0)$表示二階導(dǎo)數(shù)矩陣。若Hessian矩陣是正定的，函數(shù)的局部最小值可以通過(guò)使上面的二次型的一階導(dǎo)數(shù)等于0來(lái)獲取，我們有： [ \mathbf{x}_{\mathrm{opt}} = \mathbf{x}_0 - \mathbf{H}^{-1}\nabla f ]

這里可使用共軛梯度近似Hessian矩陣的逆矩陣。下面給出Rosenbrock函數(shù)的Hessian矩陣元素通式：

[ \begin{split} H{i,j} = \frac{\partial^2 f}{\partial x_i \partial x_j} &= 200(\delta{i,j} - 2x{i-1}\delta{i-1,j}) - 400xi(\delta{i+1,j} - 2xi\delta{i,j}) - 400\delta{i,j}(x{i+1} - xi^2) + 2\delta{i,j}, \ &= (202 + 1200xi^2 - 400x{i+1}) \delta{i,j} - 400x_i\delta{i+1,j} - 400x{i-1}\delta{i-1,j} \end{split}] 其中$i,j \in [1, N-2]$。其他邊界上的元素通式為： [ \begin{split} \frac{\partial^2 f}{\partial x_0^2} &= 1200x_0^2 - 400x_1 + 2, \ \frac{\partial^2 f}{\partial x_0 \partial x_1} = \frac{\partial^2 f}{\partial x_1 \partial x_0} &= -400x_0, \ \frac{\partial^2 f}{\partial x{N-1} \partial x{N-2}} = \frac{\partial^2 f}{\partial x{N-2} \partial x{N-1}} &= -400x_{N-2}, \ \frac{\partial^2 f}{\partial x_{N-1}^2} &= 200. \end{split}]

例如，當(dāng)$N=5$時(shí)的Hessian矩陣為：

[ \mathbf{H} =
[1200x20?400x1+2?400x0000 ?400x0202+1200x21?400x2?400x100 0?400x1202+1200x22?400x3?400x20 00?400x2202+1200x23?400x4?400x3 000?400x3200]

]

為使用牛頓共軛梯度法，我們需要提供一個(gè)計(jì)算Hessian矩陣的函數(shù)：
In [18]:

def rosen_hess(x):
    x = np.asarray(x)
    H = np.diag(-400*x[:-1],1) - np.diag(400*x[:-1],-1)
    diagonal = np.zeros_like(x)
    diagonal[0] = 1200*x[0]**2-400*x[1]+2
    diagonal[-1] = 200
    diagonal[1:-1] = 202 + 1200*x[1:-1]**2 - 400*x[2:]
    H = H + np.diag(diagonal)
    return H

In [19]:

res = opt.minimize(rosen, x_0, method='Newton-CG', jac=rosen_der, hess=rosen_hess, options={'xtol': 1e-8, 'disp': True})
print "Result of minimizing Rosenbrock function via Newton-Conjugate-Gradient algorithm (Hessian):"
print res

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 20
         Function evaluations: 22
         Gradient evaluations: 41
         Hessian evaluations: 20
Result of minimizing Rosenbrock function via Newton-Conjugate-Gradient algorithm (Hessian):
  status: 0
 success: True
    njev: 41
    nfev: 22
     fun: 1.47606641102778e-19
       x: array([ 1.,  1.,  1.,  1.,  1.])
 message: 'Optimization terminated successfully.'
    nhev: 20
     jac: array([ -3.62847530e-11,   2.68148992e-09,   1.16637362e-08,
         4.81693414e-08,  -2.76999090e-08])

對(duì)于一些大型的優(yōu)化問(wèn)題，Hessian矩陣將異常大，牛頓共軛梯度法用到的僅是Hessian矩陣和一個(gè)任意向量的乘積，為此，用戶可以提供兩個(gè)向量，一個(gè)是Hessian矩陣和一個(gè)任意向量$\mathbf{p}$的乘積，另一個(gè)是向量$\mathbf{p}$，這就減少了存儲(chǔ)的開(kāi)銷。記向量$\mathbf{p} = (p_1, \ldots, p_{N-1})$，可有

[ \mathbf{H(x)p} = \begin{bmatrix} (1200x0^2 - 400x_1 + 2)p_0 -400x_0p_1 \ \vdots \ -400x{i-1}p{i-1} + (202 + 1200x_i^2 - 400x{i+1})pi - 400x_ip{i+1} \ \vdots \ -400x{N-2}p{N-2} + 200p_{N-1} \end{bmatrix} ]

我們定義如下函數(shù)并使用牛頓共軛梯度方法尋優(yōu)：
In [20]:

def rosen_hess_p(x, p):
    x = np.asarray(x)
    Hp = np.zeros_like(x)
    Hp[0] = (1200*x[0]**2 - 400*x[1] + 2)*p[0] - 400*x[0]*p[1]
    Hp[1:-1] = -400*x[:-2]*p[:-2]+(202+1200*x[1:-1]**2-400*x[2:])*p[1:-1] \
               -400*x[1:-1]*p[2:]
    Hp[-1] = -400*x[-2]*p[-2] + 200*p[-1]
    return Hp

res = opt.minimize(rosen, x_0, method='Newton-CG', jac=rosen_der, hessp=rosen_hess_p, options={'xtol': 1e-8, 'disp': True})
print "Result of minimizing Rosenbrock function via Newton-Conjugate-Gradient algorithm (Hessian times arbitrary vector):"
print res

Optimization terminated successfully.
         Current function value: 0.000000
         Iterations: 20
         Function evaluations: 22
         Gradient evaluations: 41
         Hessian evaluations: 58
Result of minimizing Rosenbrock function via Newton-Conjugate-Gradient algorithm (Hessian times arbitrary vector):
  status: 0