共轭方向法和共轭梯度法(二次型)
文章目录
共轭方向法
共轭向量法是一种效率介于最速下降法和牛顿法之间的算法，具有以下特性：
- 对于n维的二次型函数，能够在n步之内收敛
- 不需要计算黑塞矩阵
- 不需要对黑塞矩阵求逆
共轭向量的定义
-
正交向量：
e i T e j = 0 e^T_ie_j = 0 eiTej=0 -
共轭向量
w T A v = 0 w^TAv = 0 wTAv=0
这里的A是正定对称矩阵 A = A T A = A^T A=AT记作** < w , v > A <w,v>_A <w,v>A**
本征向量相互共轭
- 有一组本征向量 x 1 − x n x_1-x_n x1−xn
- 有以下性质： A x i = λ x i Ax_i = \lambda x_i Axi=λxi
- 并且本征向量相互正交： x i T x i = 0 x_i^Tx_i = 0 xiTxi=0
- 所以有： x j T A x i = λ i x j T x i = 0 x^T_jAx_i = \lambda_i x^T_jx_i = 0 xjTAxi=λixjTxi=0
一组相互共轭向量相互独立
- 存在正定对称矩阵A，一组共轭向量： u 1 − u j u_1-u_j u1−uj
- 使用反证法，如果假设不成立，即一组共轭向量之间有线性关系，那么一定存在
∑ i = 1 n a i u i = 0 \sum_{i=1}^na_iu_i = 0 i=1∑naiui=0
其中 a i a_i ai是 常数，且不全为0
-
左乘A矩阵
A ∑ i = 1 n a i u i = ∑ i = 1 n a i A u i A\sum_{i=1}^na_iu_i = \sum_{i=1}^na_iAu_i Ai=1∑naiui=i=1∑naiAui -
左乘另外一个共轭向量 u k u_k uk
u i T ∑ i = 1 n a i A u i = ∑ i = 1 n a i u i T A u i u_i^T\sum_{i=1}^na_iAu_i = \sum_{i=1}^na_iu_i^TAu_i uiTi=1∑naiAui=i=1∑naiuiTAui -
因为 u i 和 u K u_i和u_K ui和uK是共轭向量，所以原式就变为
a k < u k , u k > = 0 a_k<u_k,u_k> =0 ak<uk,uk>=0 -
所以只能ak等于0，且k=1,2,3…n，也就是说a全为0，与原假设矛盾，所以共轭函数线性独立
一组共轭向量作为基
- 将一组向量作为基底，x可以用一组共轭向量u表示：
x = ∑ i = 1 n a i u i x = \sum_{i=1}^n a_iu_i x=i=1∑naiui
- 左乘 u k T A u^T_kA ukTA可以得到
u k T A x = ∑ i = 1 n a i u k T A u i u^T_kAx = \sum^n_{i=1}a_iu^T_kAu_i ukTAx=i=1∑naiukTAui
- 由于只有当 i ≠ k i\not=k i=k的时候， u k T A u i ≠ 0 u^T_kAu_i\not=0 ukTAui=0,所以等式可化为：
u k T A x = a k u k T A u k u^T_kAx = a_ku^T_kAu_k ukTAx=akukTAuk
a k = u k T A x u k T A u k a_k = \frac{u^T_kAx}{u_k^TAu_k} ak=ukTAukukTAx
-
将3得到的式子带入到1中可得：
x = ∑ i = 1 n u i T A x u i T A u i u i x = \sum_{i=1}^n \frac{u^T_iAx}{u_i^TAu_i}u_i x=i=1∑nuiTAuiuiTAxui线性方程组的求解
如果将上述结果用来求解线性方程组会发生什么呢？
假设Ax = b成立
x = ∑ i = 1 n u i T A x u i T A u i u i ⟹ x = ∑ i = 1 n u i T b u i T A u i u i x = \sum_{i=1}^n \frac{u^T_iAx}{u_i^TAu_i}u_i \Longrightarrow x = \sum_{i=1}^n \frac{u^T_ib}{u_i^TAu_i}u_i x=i=1∑nuiTAuiuiTAxui⟹x=i=1∑nuiTAuiuiTbui应用到牛顿法中
之前牛顿法的迭代式 X k = X k − 1 − H − 1 g X_k = X_{k-1}-H^{-1}g Xk=Xk−1−H−1g其中g是梯度。
牛顿法的一个大问题就是H求逆的工作量太大，耗费太多资源，我们怎么能不求逆的情况下迭代牛顿法呢？这就是我们要解决的问题
我们令 d = − H − 1 g ⟹ H d = − g d = -H^{-1}g\Longrightarrow Hd = -g d=−H−1g⟹Hd=−g
根据上一步的推导，我们可以得到： d = ∑ i = 1 n − u i T g u i T H u i u i d = \sum_{i=1}^n -\frac{u^T_ig}{u_i^THu_i}u_i d=∑i=1n−uiTHuiuiTgui,
式子分母是一个标量，分母也是一个标量，这样我们就巧妙的规避掉复杂的求逆过程。且H也是一个正定并且对称的矩阵
共轭方向法在二次型函数的应用
二次型目标函数可以写作：
f
(
x
)
=
1
2
X
T
H
X
−
b
T
X
f(x) = \frac12X^THX-b^TX
f(x)=21XTHX−bTX
梯度（一阶导）可以写作
g
k
=
H
X
K
−
b
g_k = HX_K-b
gk=HXK−b
共轭方向法
X
∗
=
X
1
−
H
−
1
g
(
X
1
)
⟹
X
∗
=
X
1
−
∑
i
=
1
n
u
i
T
g
(
X
1
)
u
i
T
H
u
i
u
i
为
X^* = X_1-H^{-1}g(X_1) \Longrightarrow X^* = X_1-\sum_{i=1}^n \frac{u^T_ig(X_1)}{u_i^THu_i}u_i为
X∗=X1−H−1g(X1)⟹X∗=X1−i=1∑nuiTHuiuiTg(X1)ui为
这里迭代的方向d为 u i u_i ui(和上面的d没有关系),方向d根据假设，应该是一组互相共轭的方向向量
步长 a i = d i T g ( X 1 ) d i T H d i a_i = \frac{d^T_ig(X_1)}{d_i^THd_i} ai=diTHdidiTg(X1)
将共轭方向法写成一般形式就是：
X
∗
=
X
0
+
a
1
d
1
+
a
2
d
2
+
a
3
d
3
.
.
.
+
a
n
d
n
X^* = X_0+a_1d_1+a_2d_2+a_3d_3...+a_nd_n
X∗=X0+a1d1+a2d2+a3d3...+andn
X
∗
−
X
0
=
a
1
d
1
+
a
2
d
2
+
a
3
d
3
.
.
.
+
a
n
d
n
d
k
T
H
(
X
∗
−
X
0
)
=
a
1
d
k
T
H
d
1
+
a
2
d
k
T
H
d
2
.
.
.
+
a
n
d
k
T
H
d
n
d
k
T
H
(
X
∗
−
X
0
)
=
a
k
d
k
T
H
d
k
a
k
=
d
k
T
H
(
X
∗
−
X
0
)
d
k
T
H
d
k
X
∗
−
X
0
=
(
X
∗
−
X
k
)
+
(
X
k
−
X
1
)
a
k
=
d
k
T
H
[
(
X
∗
−
X
k
)
+
(
X
k
−
X
1
)
]
d
k
T
H
d
k
a
k
=
−
d
k
T
g
k
d
k
T
H
d
k
X^*-X_0 = a_1d_1+a_2d_2+a_3d_3...+a_nd_n\\ d_k^TH(X^*-X_0) = a_1d_k^THd_1+a_2d_k^THd_2...+a_nd_k^THd_n\\ d_k^TH(X^*-X_0) = a_kd_k^THd_k\\ a_k = \frac{d_k^TH(X^*-X_0)}{d_k^THd_k}\\ X^*-X_0 = (X^* - X_k) + (X_k - X_1)\\ a_k = \frac{d_k^TH[(X^* - X_k) + (X_k - X_1)]}{d_k^THd_k}\\ a_k = -\frac{d^T_kg_k}{d_k^THd_k}\\
X∗−X0=a1d1+a2d2+a3d3...+andndkTH(X∗−X0)=a1dkTHd1+a2dkTHd2...+andkTHdndkTH(X∗−X0)=akdkTHdkak=dkTHdkdkTH(X∗−X0)X∗−X0=(X∗−Xk)+(Xk−X1)ak=dkTHdkdkTH[(X∗−Xk)+(Xk−X1)]ak=−dkTHdkdkTgk
这里算出来的
a
k
=
−
d
k
T
g
k
d
k
T
H
d
k
a_k = -\frac{d^T_kg_k}{d_k^THd_k}
ak=−dkTHdkdkTgk与上面的
a
i
=
d
i
T
g
(
X
0
)
d
i
T
H
d
i
a_i = \frac{d^T_ig(X_0)}{d_i^THd_i}
ai=diTHdidiTg(X0)其实是等价的
证明：
d
k
T
g
k
=
d
k
T
(
H
X
k
−
b
)
=
d
k
T
[
H
(
X
0
+
a
1
d
1
+
a
2
d
2
+
a
3
d
3
.
.
.
+
a
n
d
n
)
−
b
]
=
d
k
T
(
H
X
0
−
b
)
=
d
k
T
g
0
d^T_kg_k = d^T_k(HX_k -b)\\ = d^T_k[H(X_0+a_1d_1+a_2d_2+a_3d_3...+a_nd_n)-b]\\ = d^T_k(HX_0 -b)\\ =d^T_kg_0
dkTgk=dkT(HXk−b)=dkT[H(X0+a1d1+a2d2+a3d3...+andn)−b]=dkT(HX0−b)=dkTg0
所以这里不论是使用
a
k
=
−
d
k
T
g
k
d
k
T
H
d
k
a_k = -\frac{d^T_kg_k}{d_k^THd_k}
ak=−dkTHdkdkTgk还是
a
i
=
d
i
T
g
(
X
0
)
d
i
T
H
d
i
a_i = \frac{d^T_ig(X_0)}{d_i^THd_i}
ai=diTHdidiTg(X0)都是可以的
共轭方向法的定理
-
对于任意的初始点，共轭方向法都能在n次迭代之内收敛到唯一的全局极小点。
-
在算法中，都有
g ( k + 1 ) T d i = 0 , i = 0 , 1..... k g(k+1)^Td^i=0,i=0,1.....k g(k+1)Tdi=0,i=0,1.....k该定理被称为扩张子空间定理，我一直不理解一个向量怎么同时垂直不在一个平面上的向量，能请有看到的大佬解释一下吗？？？？？（博主的数学是真的不行）
共轭方向法的效率很高，但是必须能够给定一组共轭的向量，接下来我们就将要解决这个问题。
二次型共轭梯度法
Gram-Schmidt过程
-
构建正交基
假如存在一组基底v（不一定正交),那可以通过以下方法构造一组正交基
u 1 = v 1 u_1 = v_1 u1=v1u 2 = v 2 − < v 2 , u 1 > < u 1 , u 1 > u 1 u_2 = v_2 - \frac{<v_2,u_1>}{<u_1,u_1>}u_1 u2=v2−<u1,u1><v2,u1>u1
$$
u_3 = v_3 - \frac{<v_3,u_1>}{<u_1,u_1>}u_1- \frac{<v_3,u_2>}{<u_2,u_2>}u_2
$$
u k = v k − ∑ i = 1 k − 1 β i k u i u_k = v_k - \sum_{i=1}^{k-1} \beta^k_iu_i uk=vk−i=1∑k−1βikui
其中 β i k = < v k , u i > < u i , u i > \beta^k_i = \frac{<v_k,u_i>}{<u_i,u_i>} βik=<ui,ui><vk,ui>
通过数学归纳法可以证明以上向量组u两两正交
- \frac{<v_3,u_2>}{<u_2,u_2>}u_2
-
构建共轭基
u k = v k − ∑ i = 1 k − 1 β i k u i u_k = v_k - \sum_{i=1}^{k-1} \beta^k_iu_i uk=vk−i=1∑k−1βikui
其中 β i k = < v k , u i > H < u i , u i > H \beta^k_i= \frac{<v_k,u_i>_H}{<u_i,u_i>_H} βik=<ui,ui>H<vk,ui>H -
共轭梯度法
-
共轭梯度法：初始搜索方向向量设置为-g，梯度负方向
-
一阶导条件： g = H X ∗ − b = 0 g = HX^*-b=0 g=HX∗−b=0 ，令 r k = b − H X K = − g k r_k = b - HX_K = -g_k rk=b−HXK=−gk用来判断梯度是否足够小，来终止迭代
-
那么搜索方向对迭代点的梯度进行构建，使其变成一组相互共轭的方向向量组：
d k = − g k − ∑ i = 1 k − 1 β i k d i d_k = -g_k - \sum_{i=1}^{k-1}\beta_i^kd_i dk=−gk−∑i=1k−1βikdi,其中 β i k = − < g k , d i > H < d i , d i > H \beta^k_i= -\frac{<g_k,d_i>_H}{<d_i,d_i>_H} βik=−<di,di>H<gk,di>H
-
d k = r k − ∑ i = 1 k − 1 β i k d i d_k = r_k - \sum_{i=1}^{k-1}\beta_i^kd_i dk=rk−∑i=1k−1βikdi,其中 β i k = < r k , d i > H < d i , d i > H \beta^k_i= \frac{<r_k,d_i>_H}{<d_i,d_i>_H} βik=<di,di>H<rk,di>H
-
搜索步长： a k = d K T r K d k T H d k a_k = \frac{d^T_Kr_K}{d^T_kHd_k} ak=dkTHdkdKTrK
-
残差e： e k = x k − x ∗ e_k = x_k - x^* ek=xk−x∗ , x ∗ x^* x∗是极小值，可得残差的性质：
H e k = H x k − H x ∗ = H x k − b = g k = − r k He_k = Hx_k - Hx^*\\ = Hx_k -b\\ =g_k=-r_k Hek=Hxk−Hx∗=Hxk−b=gk=−rk -
残差与搜索方向共轭：
x k = x 1 + a 1 d 1 + . . . a k − 1 d k − 1 x ∗ = x 1 + a 1 d 1 + . . . a n d n e = ∑ i = k n a i d i d i T H e k = 0 , i = 1 , 2 , . . . . k − 1 d i T g k = 0 , i = 1 , 2 , . . . . k − 1 x_k= x_1 + a_1d_1+ ... a_{k-1}d_{k-1}\\ x^* = x_1 + a_1d_1+ ... a_nd_n\\ e = \sum_{i=k}^n a_id_i\\ d_i^THe_k = 0,i=1,2,....k-1\\ d_i^Tg_k = 0,i=1,2,....k-1 xk=x1+a1d1+...ak−1dk−1x∗=x1+a1d1+...andne=i=k∑naididiTHek=0,i=1,2,....k−1diTgk=0,i=1,2,....k−1
也就是说残差与之前所有的搜索方向共轭，并且当前梯度与之前所有搜索方向正交
-
梯度方向相互正交
d i T g k = 0 , i = 1 , 2... k d_i^Tg_k = 0,i=1,2...k diTgk=0,i=1,2...k
梯度方向 d i = r i − ∑ = 1 k − 1 β j k d j d_i = r_i - \sum_{=1}^{k-1}\beta_j^kd_j di=ri−∑=1k−1βjkdj
将梯度方向代入可以得到：
g k T d i = g k T r i − ∑ j = 1 k − 1 β j g k T d j = 0 g^T_kd_i = g_k^Tr_i - \sum_{j=1}^{k-1}\beta_jg^T_kd_j = 0 gkTdi=gkTri−j=1∑k−1βjgkTdj=0
可以得到：
g k T r i = 0 = − g k T g i g_k^Tr_i = 0 = -g_k^Tg_i gkTri=0=−gkTgi -
当前点的方向与之前所有点的梯度相互共轭：
x k + 1 = x k + a k d k H x k + 1 = H x k + a k H d k b − H x k + 1 = b − H x k − a k H d k r k + 1 = r k − a k H d k x_{k+1} = x_k +a_kd_k\\ Hx_{k+1} = Hx_k + a_kHd_k\\ b -Hx_{k+1} = b-Hx_k - a_kHd_k\\ r_{k+1} = r_k-a_kHd_k xk+1=xk+akdkHxk+1=Hxk+akHdkb−Hxk+1=b−Hxk−akHdkrk+1=rk−akHdk
这样就得到了一个 r k + 1 r_{k+1} rk+1的迭代式
H d k = 1 a k ( r k − r k + 1 ) r i T H d k = 1 a k ( r i T r k − r i T r k + 1 ) Hd_k = \frac1{a_k}(r_k-r_{k+1})\\ r^T_iHd_k = \frac1{a_k}(r^T_ir_k-r^T_ir_{k+1})\\ Hdk=ak1(rk−rk+1)riTHdk=ak1(riTrk−riTrk+1)r i T H d k = { 1 a k r k T r k i = k − 1 a k r k + 1 T r k + 1 i = k + 1 0 o t h e r w i s e r^T_iHd_k = \begin{cases} \frac1{a_k}r^T_kr_k & i=k\\ -\frac1{a_k}r^T_{k+1}r_{k+1} & i = k+1\\ 0 & otherwise \end{cases} riTHdk=⎩⎪⎨⎪⎧ak1rkTrk−ak1rk+1Trk+10i=ki=k+1otherwise
这个结果也就是说明，当前点的方向与除了当前点和下一个点以外的所有点的梯度相互共轭，这边与上面第五点的结论有区别，5的结论是当前迭代点的梯度与之前所有的搜索方向正交
-
改进Gram-Schmidt过程
Gram-Schmidt过程：
d k = r k − ∑ i = 1 k − 1 β i k d i d_k = r_k - \sum_{i=1}^{k-1} \beta^k_id_i dk=rk−i=1∑k−1βikdi
其中 β i k = − < g k , d i > H < d i , d i > H \beta^k_i= -\frac{<g_k,d_i>_H}{<d_i,d_i>_H} βik=−<di,di>H<gk,di>H如果我们直接使用的话，需要记住之前计算的所有搜索方向,会浪费很多资源。所以我们对Gram-Schmidt过程进行改进
β i k = − < g k , d i > H < d i , d i > H 因 为 只 有 i = k − 1 的 时 候 β ≠ 0 根 据 7 的 结 论 d k = r k − β k − 1 d k − 1 其 中 β k − 1 = r k T H d k − 1 d k − 1 T H d k − 1 \beta^k_i= -\frac{<g_k,d_i>_H}{<d_i,d_i>_H}\\ 因为只有i=k-1的时候\beta\not=0\\ 根据7的结论 \\ d_k = r_k - \beta_{k-1}d_{k-1}\\ 其中\beta_{k-1} = \frac{r_k^THd_{k-1}}{d^T_{k-1}Hd_{k-1}} βik=−<di,di>H<gk,di>H因为只有i=k−1的时候β=0根据7的结论dk=rk−βk−1dk−1其中βk−1=dk−1THdk−1rkTHdk−1
这样在我们就不需要存储之前计算的方向向量，节约了资源。 -
迭代过程
-
初始条件：
g 1 = H x 1 − b d 1 = r 1 = − g 1 g_1 = Hx_1-b\\ d_1 = r_1 = -g_1 g1=Hx1−bd1=r1=−g1 -
迭代关系
a k = d K T r k d k T H d k x k + 1 = x k + a k d k r k + 1 = r k − a k H d k β k = r k + 1 T H d k d k T H d k d k + 1 = r k + 1 − β k d k a_k = \frac{d^T_Kr_k}{d^T_kHd_k} \\x_{k+1} = x_k + a_kd_k \\r_{k+1} = r_k - a_kHd_k \\ \beta_{k} = \frac{r_{k+1}^THd_{k}}{d^T_{k}Hd_{k}}\\d_{k+1} = r_{k+1} - \beta_kd_k ak=dkTHdkdKTrkxk+1=xk+akdkrk+1=rk−akHdkβk=dkTHdkrk+1THdkdk+1=rk+1−βkdk
-
-
二次型共轭梯度法代码实现（Julia）
using LinearAlgebra
using Gadfly
println("----------------程序开始-----------------")
function ConjGradQuad(b,
H,
start;
accuracy = 0.01,
debug = false,
plot = false)
## 判断输入是否合法
if ndims(H) !=2
error("黑塞矩阵应该是个矩阵")
end
nr,nc = size(H)
if nr != nc
error("黑塞矩阵的行数应该等于列数")
end
## 将H强制转化为对称矩阵（按照上三角矩阵来复制下三角矩阵）
h = Symmetric(H)
if !(h == H)
println("warning:输入的黑塞矩阵不是对称的，已自动转化为对称矩阵 ")
end
## 为下面的作图做准备
point = []
push!(point,start)
x1 = start
g1 = h * x1 - b
d1 = r1 = -1 .* g1
for i in 1:nc
α = d1' * r1 /(d1' *h *d1)
xk = x1 + α .*d1
rk = r1 = α .* h *d1
## 输出过程
if debug
println("x:",xk,"\tr:",rk)
end
## 如果要画图的话
if plot
push!(point,xk)
end
## 用梯度足够小（一阶导条件）来结束煦暖
if norm(rk) < accuracy
return xk,rk,point,i
end
βk = rk' *h *d1 / (d1' *h *d1)
dk = rk - βk .* d1
## upload
d1 = dk
r1 = rk
x1 = xk
end
println("最大迭代数：",nc,"，无法达到要求的梯度精度")
return x1,r1,point
end
## 测试
## 反悔了极小点，梯度，绘画的array以及迭代的次数
x1,r1,point = ConjGradQuad([4,4],
[5 -3;-3 5],
[8,6],
plot = true,
debug = false)
println("----------------最终结果-----------------")
println("x:",x1)
println("r:",r1)
println(point)
## 一定要记得使用小数
low = -0.0 ##下线
up = 20.0
## 点数等于迭代次数+1
n = length(point)
rosen = layer((x,y) -> 5/2*(x^2 + y^2) - 3*x*y -4*(x+y),low,up,low,up,Geom.contour())
steps = layer(
## 最后一点不作为起点
x = convert(Vector{Float64},[point[i][1] for i in 1:n-1]),
y = convert(Vector{Float64},[point[i][2] for i in 1:n-1]),
## 第一点不作为终点
xend = convert(Vector{Float64},[point[i][1] for i in 2:n]),
yend = convert(Vector{Float64},[point[i][2] for i in 2:n]),
Geom.point,Geom.vector(filled = true))
plot(rosen,steps,Scale.x_continuous(minvalue = low,maxvalue = up),
Scale.y_continuous(maxvalue=up,minvalue=low))
引用包
using LinearAlgebra
using Gadfly
构造函数
需要的参数有b、H、start初始点、默认精度是0.01,。
function ConjGradQuad(b,
H,
start;
accuracy = 0.01,
debug = false,
plot = false)
end
判断输入的函数是否合法
## 判断输入是否合法
if ndims(H) !=2
error("黑塞矩阵应该是个矩阵")
end
nr,nc = size(H)
if nr != nc
error("黑塞矩阵的行数应该等于列数")
end
## 将H强制转化为对称矩阵（按照上三角矩阵来复制下三角矩阵）
h = Symmetric(H)
if !(h == H)
println("warning:输入的黑塞矩阵不是对称的，已自动转化为对称矩阵 ")
end
初始准备
## 为下面的作图做准备
point = []
push!(point,start)
x1 = start
g1 = h * x1 - b
d1 = r1 = -1 .* g1
循环迭代
for i in 1:nc
α = d1' * r1 /(d1' *h *d1)
xk = x1 + α .*d1
rk = r1 = α .* h *d1
## 用梯度足够小（一阶导条件）来结束循环
if norm(rk) < accuracy
return xk,rk,point,i
end
βk = rk' *h *d1 / (d1' *h *d1)
dk = rk - βk .* d1
## upload
d1 = dk
r1 = rk
x1 = xk
end
迭代画图需要使用到的array
point是一个用来存放迭代点的array，在数据准备的时候已经将初始迭代点添加进point，只有需要绘图，也就是plot参数为真时，才会迭代这个array
if plot
push!(point,xk)
end
构造测试用例
## 测试
## 返回了极小点，梯度，绘画的array
x1,r1,point = ConjGradQuad([4,4],
[5 -3;-3 5],
[8,6],
plot = true,
debug = false)
输出迭代结果
println("----------------最终结果-----------------")
println("x:",x1)
println("r:",r1)
将迭代的步骤画出来
这里需要注意一个点，那就是如果layer函数里面一旦有一个值是小数，那么所有的参数就一定要是浮点型的，否则会报错。（我也不知道为什么）所以为了方便起见，最好就全部都是小数。
为了让每一个参数都是浮点型，使用了**convert()**函数进行数值转换
## 画出point看看效果
println(point)
## 一定要记得使用小数
low = -0.0 ##下限
up = 20.0 ## 上限
n = length(point)
rosen = layer((x,y) -> 5/2*(x^2 + y^2) - 3*x*y -4*(x+y),low,up,low,up,Geom.contour())
steps = layer(
## convert()函数将数据转化成了浮点型，这样子才不会报错
## 最后一点不作为起点
x = convert(Vector{Float64},[point[i][1] for i in 1:n-1]),
y = convert(Vector{Float64},[point[i][2] for i in 1:n-1]),
## 第一点不作为终点
xend = convert(Vector{Float64},[point[i][1] for i in 2:n]),
yend = convert(Vector{Float64},[point[i][2] for i in 2:n]),
Geom.point,Geom.vector(filled = true))
plot(rosen,steps,Scale.x_continuous(minvalue = low,maxvalue = up),
Scale.y_continuous(maxvalue=up,minvalue=low))
这两个函数原本是使用Gadfly包的plot()函数进行绘制，但是为了可以把目标函数的图像和迭代点的图像合并起来，所以改成了使用**layer()**代替ploy()将绘画的信息存储起来。然后再统一绘画出来。
rosen = layer((x,y) -> 5/2*(x^2 + y^2) - 3*x*y -4*(x+y),low,up,low,up,Geom.contour())
steps = layer(
## convert()函数将数据转化成了浮点型，这样子才不会报错
## 最后一点不作为起点
x = convert(Vector{Float64},[point[i][1] for i in 1:n-1]),
y = convert(Vector{Float64},[point[i][2] for i in 1:n-1]),
## 第一点不作为终点
xend = convert(Vector{Float64},[point[i][1] for i in 2:n]),
yend = convert(Vector{Float64},[point[i][2] for i in 2:n]),
Geom.point,Geom.vector(filled = true))
- Geom.point表示将每个点画出来
- **Geom.vector()**表示画出箭头，filled = true表示画的是实心箭头
将目标函数和迭代步画出来
plot(rosen,
steps,
Scale.x_continuous(minvalue = low,maxvalue = up),
Scale.y_continuous(maxvalue=up,minvalue=low))
绘画结果
1219
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
