1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 1. 目录2. 模拟方法与误差 2.1. 模拟方法 2.2. 误差分析 2.2.1. 向前误差和向后误差 Forward Error and Backward Error 2.2.2. 敏感性与病态性 Sensitivity and Conditioning 2.2.3. 条件数 Conditional Number 2.2.3.1. 向量的条件数 2.2.3.1. 方阵的条件数 2.2.3.2. 误差限 2.2.4. 截断误差和舍入误差 Truncation Error and Round-off Error 3. 线性方程组的求解 3.1. 矩阵的秩 3.2. 方程组解的个数 3.2.1. 非奇异矩阵 Nonsingular Matrix 3.3. 范数/模 Norm 3.3.1. 向量的范数 3.3.1.1. Euclidean Norm 3.3.1.2. H-Norm 3.3.2. 矩阵的范数 3.3.2.1. p-范数诱导的矩阵范数 3.3.2.2. Frobenius Norm 3.4. 求解线性方程组的方法 3.5. Gauss 消去法 3.5.1. 初等消去矩阵 3.5.2. Gauss 消去法内容 3.6. 共轭梯度法 Conjugate Gradient Method 3.6.1. 先置知识 3.6.2. 矩阵A共轭向量 3.6.3. 当已知某个矩阵A共轭向量集合时 3.6.4. 共轭梯度法的收敛性 3.6.5. 迭代系数alpha的设计 3.6.6. 如何产生矩阵A共轭向量集合{p_i} 3.6.6.1. Expanding Subspace Minimization 3.6.6.2. 基于矩阵A共轭向量集合{p_i}特性产生{p_i} 3.6.7. 共轭梯度法算法步骤 3.7. 迭代法 Iterative Method (待更新) 3.7.1. 雅可比迭代法 Jacobi Method (待更新) 3.7.2. 高斯-塞德尔迭代法 Gauss-Seidel Method (待更新) 4. 线性最小二乘问题 (待更新)5. 无约束优化 5.1. 优化基础 5.2. 下降方法的收敛性 5.2.1. 关于收敛 5.2.2. 部分下降方法的收敛性 5.3. 最速下降法 Steepest Descent Method 5.4. 牛顿法 Newton-Raphson Method 5.4.1. 收敛性 5.5. 实用牛顿法 Practical Newton Method 5.5.1. 信赖域方法 Trust-Region Newton Methods 5.6. 拟牛顿法 Quasi-Newton Method 5.6.1. 拟牛顿条件/割线条件 5.6.2. DFP方法 5.6.2.1. DFP算法步骤 5.6.3. BFGS方法 5.6.3.1. BFGS算法步骤 6. 关于学习中的一些补充 (待更新)
计算模拟重要方法：迭代法、逼近法
一般策略：
用有限维空间代替无限的空间。
用有限过程代替无限过程。
用代数方程代替微分方程。
用线性问题代替非线性问题。
用低阶方程组代替高阶方程组。
用简单函数代替复杂函数。
用简单结构的矩阵代替一般矩阵。
近似的来源：
建模：这个过程可能简化，或者和忽略问题或者系统的某些物理特性。
经验测量：实验设备精度所限。
前面的计算：输入的数据可能是中间数据，其结果只能是近似值。
绝对误差 = 近似值 - 真值
相对误差 = 绝对误差 / 真值
总误差 = 计算误差 + 数据传播误差 ，比如
Δ f ( x ) = f ^ ( x ^ ) − f ( x ) = ( f ^ ( x ^ ) − f ( x ^ ) ) ⏟ C o m p u t a t i o n a l E r r o r + ( f ( x ^ ) − f ( x ) ) ⏟ P r o p a g a t e d D a t a E r r o r \Delta f(x)=\hat f(\hat x)-f(x)=\underbrace{(\hat f(\hat x)-f(\hat x))}_{ComputationalError}+\underbrace{(f(\hat x)-f(x))}_{PropagatedDataError}
Δ f ( x ) = f ^ ( x ^ ) − f ( x ) = C o m p u t a t i o n a l E r r o r ( f ^ ( x ^ ) − f ( x ^ ) ) + P r o p a g a t e d D a t a E r r o r ( f ( x ^ ) − f ( x ) ) 若有 y = f ( x ) a n d y ^ = f ( x ^ ) y=f(x)\quad and \quad\hat y=f(\hat x) y = f ( x ) a n d y ^ = f ( x ^ ) Δ y = y ^ − y \Delta y=\hat y-y Δ y = y ^ − y Δ x = x ^ − x \Delta x=\hat x-x Δ x = x ^ − x
fig1.2
f ( x ^ ) = f ^ ( x ) f(\hat x)=\hat f(x) f ( x ^ ) = f ^ ( x ) x ^ \hat x x ^
Note that the equality $\hat f(x)=f(\hat x)$ is due to the choice of $\hat x$; indeed, this requirement *defines* $\hat x$.
在某些情况下，即便计算是完全精确的，问题的解对数据的扰动 也可能是高度敏感的。
The qualitative notion of *insensitive*, and its quantitative measure, called *conditioning*, are concerned with propagated data error, i.e., the effects on the solution of perturbations in the input data.
如果输入数据的相对变化对解的相对变化影响相对适中，则是不敏感 (insensitive ) 的，或者良态 (well-conditioned ) 的。
如果影响很大，则是敏感 (sensitive ) 的，或者病态 (ill-conditioned ) 的。
解的相对变化与输入数据的相对变化的比值称为条件数。(一个问题的条件数远远大于1，则是病态的。)
C o n d i t i o n n u m b e r = ∣ ( f ( x ^ ) − f ( x ) ) / f ( x ) ∣ ∣ ( x ^ − x ) / x ∣ = ∣ y ^ − y y ∣ ∣ x ^ − x x ∣ = ∣ Δ y / y Δ x / x ∣ Condition\quad number=\frac{|(f(\hat x)-f(x))/f(x)|}{|(\hat x-x)/x|}=\frac{|\frac{\hat y-y}{y}|}{|\frac{\hat x-x}{x}|}=|\frac{\Delta y/y}{\Delta x/x}|
C o n d i t i o n n u m b e r = ∣ ( x ^ − x ) / x ∣ ∣ ( f ( x ^ ) − f ( x ) ) / f ( x ) ∣ = ∣ x x ^ − x ∣ ∣ y y ^ − y ∣ = ∣ Δ x / x Δ y / y ∣
∣ Δ y / y ∣ ⏟ R e l a t i v e F o r w a r d E r r o r = C o n d i t i o n n u m b e r × ∣ Δ x / x ∣ ⏟ R e l a t i v e B a c k w a r d E r r o r \underbrace{|\Delta y/y|}_{Relative\quad Forward\quad Error} = Condition\quad number \times \underbrace{|\Delta x/x|}_{Relative\quad Backward\quad Error}
R e l a t i v e F o r w a r d E r r o r ∣ Δ y / y ∣ = C o n d i t i o n n u m b e r × R e l a t i v e B a c k w a r d E r r o r ∣ Δ x / x ∣ 取$\hat f(x)-f(x)=f(\hat x)-f(x)\simeq f'(x)\Delta x$，有条件数≈$|\frac{xf'(x)}{f(x)}|$
x ⃗ = ( x 1 , x 2 , ⋯ , x n ) T , y ⃗ = F ( x ⃗ ) = ( f 1 ( x ⃗ ) , f 2 ( x ⃗ ) , ⋯ , f n ( x ⃗ ) ) T \vec x=(x_1,x_2,\cdots,x_n)^T,\qquad \vec y=\boldsymbol{F}(\vec x)=(f_1(\vec x),f_2(\vec x),\cdots,f_n(\vec x))^T
x = ( x 1 , x 2 , ⋯ , x n ) T , y = F ( x ) = ( f 1 ( x ) , f 2 ( x ) , ⋯ , f n ( x ) ) T
Δ y i = ∑ j = 1 n ∂ f i ( x ⃗ ) ∂ x j Δ x j \Delta y_i=\sum_{j=1}^n\frac{\partial f_i(\vec x)}{\partial x_j}\Delta x_j
Δ y i = j = 1 ∑ n ∂ x j ∂ f i ( x ) Δ x j 由于
ϵ y i = Δ y i f i ( x ⃗ ) , ϵ x j = Δ x j x j \epsilon_{y_i}=\frac{\Delta y_i}{f_i(\vec x)},\quad \epsilon_{x_j}=\frac{\Delta x_j}{x_j}
ϵ y i = f i ( x ) Δ y i , ϵ x j = x j Δ x j 所以
ϵ y i = 1 f i ( x ⃗ ) ∑ j = 1 n ∂ f i ( x ⃗ ) ∂ x j x j ϵ x j \epsilon_{y_i}=\frac{1}{f_i(\vec x)}\sum_{j=1}^n\frac{\partial f_i(\vec x)}{\partial x_j}x_j \epsilon_{x_j}
ϵ y i = f i ( x ) 1 j = 1 ∑ n ∂ x j ∂ f i ( x ) x j ϵ x j 即条件数
C o n d i t i o n n u m b e r = ( x 1 f 1 ( x ⃗ ) ∂ f 1 ( x ⃗ ) ∂ x 1 x 2 f 1 ( x ⃗ ) ∂ f 1 ( x ⃗ ) ∂ x 2 ⋯ x 1 f 2 ( x ⃗ ) ∂ f 2 ( x ⃗ ) ∂ x 1 x 2 f 2 ( x ⃗ ) ∂ f 2 ( x ⃗ ) ∂ x 2 ⋯ ⋮ ⋮ ⋱ ) Condition\quad number=\left(
\begin{array}{ccc}
\frac{x_1}{f_1(\vec x)}\frac{\partial f_1(\vec x)}{\partial x_1} & \frac{x_2}{f_1(\vec x)}\frac{\partial f_1(\vec x)}{\partial x_2} & \cdots\\
\frac{x_1}{f_2(\vec x)}\frac{\partial f_2(\vec x)}{\partial x_1} & \frac{x_2}{f_2(\vec x)}\frac{\partial f_2(\vec x)}{\partial x_2} & \cdots\\
\vdots & \vdots & \ddots
\end{array} \right) C o n d i t i o n n u m b e r = ⎝ ⎜ ⎜ ⎛ f 1 ( x ) x 1 ∂ x 1 ∂ f 1 ( x ) f 2 ( x ) x 1 ∂ x 1 ∂ f 2 ( x ) ⋮ f 1 ( x ) x 2 ∂ x 2 ∂ f 1 ( x ) f 2 ( x ) x 2 ∂ x 2 ∂ f 2 ( x ) ⋮ ⋯ ⋯ ⋱ ⎠ ⎟ ⎟ ⎞ 非奇异矩阵的条件数定义为c o n d ( A ) = ∥ A ∥ ⋅ ∥ A − 1 ∥ cond(\mathbf A)=\|\mathbf A\|\cdot \|\mathbf A^{-1}\| c o n d ( A ) = ∥ A ∥ ⋅ ∥ A − 1 ∥
***Proof:***
在有误差的情况下，有$\mathbf A(\vec x+\delta\vec x)=\vec b+\delta \vec b$。
因此$\|\delta \vec x\|=\|\mathbf A^{-1}\delta \vec b\|\leq \|\mathbf A^{-1}\|\|\delta \vec b\|$.
又$\|\vec b\|=\|\mathbf A\vec x\|\leq\|A\|\cdot\|x\|\Rightarrow \frac{1}{\|\vec x\|}\leq \frac{\|\mathbf A\|}{\|\vec b\|}$.
两式相乘，得到$\frac{\delta \vec x}{\|\vec x\|}\leq\|A\|\|A^{-1}\| \frac{\|\delta \vec b\|}{\|\vec b\|}$.
故条件数为$cond(\mathbf A)=\|\mathbf A\|\cdot \|\mathbf A^{-1}\|$。
方程组解的相对变化，可用条件数乘以问题中数据的相对变化来界定。δ x ⃗ ∥ x ⃗ ∥ ≤ c o n d ( A ) ⋅ ( ∥ Δ A ∥ ∥ A ∥ + ∥ δ b ⃗ ∥ ∥ b ⃗ ∥ ) \frac{\delta \vec x}{\|\vec x\|}\leq cond(A) \cdot(\frac{\|\Delta \mathbf A\|}{\|\mathbf A\|}+ \frac{\|\delta \vec b\|}{\|\vec b\|}) ∥ x ∥ δ x ≤ c o n d ( A ) ⋅ ( ∥ A ∥ ∥ Δ A ∥ + ∥ b ∥ ∥ δ b ∥ )
***Proof:***
取$\mathbf A(t)=\mathbf A+t\Delta\mathbf A,\quad\vec b(t)=\vec b+t\delta\vec b$, 并考虑线性方程组$\mathbf A(t)\vec x(t)=\vec b(t)$ 的解 $\vec x(t)=\vec x+t\delta\vec x$。求导并使$t\to 0$后取范数并整理即得该不等式。
计算误差（计算过程产生的误差）可以分为截断误差（级数的截断，迭代的中止）和舍入误差（有限精度）。
A x ⃗ = b ⃗ A ∈ R m × n , b ⃗ ∈ R m \mathbf A\vec x=\vec b\quad \mathbf A \in \mathbb R^{m\times n}, \vec b\in \mathbb R^{m} A x = b A ∈ R m × n , b ∈ R m
矩阵的 秩
是矩阵列向量生成的向量空间的维数 (Bourbaki, Algebra )
是极大线性无关组的元素向量个数
是矩阵作为线性映射变换矩阵生成的像的维度
矩阵的列秩和行秩总是相等的，因此它们可以简单地称作矩阵的秩。
对于方阵A，若满足下列等价条件之一则为非奇异；否则为奇异。
A的逆矩阵存在
行列式d e t ( A ) ≠ 0 det(\mathbf A)\neq 0 d e t ( A ) = 0
n阶矩阵的秩R a n k ( A ) = n Rank(\mathbf A)=n R a n k ( A ) = n
对任何非零向量 z ⃗ ≠ 0 , A z ⃗ ≠ 0 \vec z\neq 0,\mathbf A\vec z\neq 0 z = 0 , A z = 0
所以
若矩阵A非奇异，无论b取何值，方程组 总有唯一解
若矩阵A奇异，则解的个数由b取值决定，可能无解，也可能由无数解，且必存在非零向量z使得A z ⃗ = 0 \mathbf A\vec z=0 A z = 0
有解的充要条件 ： R a n k ( A ) = R a n k [ A , b ⃗ ] Rank(\mathbf A)=Rank[\mathbf A,\vec b] R a n k ( A ) = R a n k [ A , b ]
Norm是具有“长度”概念的函数。在线性代数、泛函分析及相关的数学领域，是一个函数，其为向量空间内的所有向量赋予非零的正长度或大小，满足
半正定性；n o r m ( v ) ≥ 0 norm(v)\geq 0 n o r m ( v ) ≥ 0
齐次性；n o r m ( a v ) = ∣ a ∣ ⋅ n o r m ( v ) norm(av)=|a|\cdot norm(v) n o r m ( a v ) = ∣ a ∣ ⋅ n o r m ( v )
三角不等式成立；n o r m ( v 1 ) + n o r m ( v 2 ) ≥ n o r m ( v 1 + v 2 ) norm(v_1)+norm(v_2)\geq norm(v_1+v_2) n o r m ( v 1 ) + n o r m ( v 2 ) ≥ n o r m ( v 1 + v 2 )
正定性（不满足则为半范数） ；n o r m ( v ) = 0 i f f v ⃗ = 0 ⃗ norm(v)= 0\quad iff\quad \vec v=\vec 0 n o r m ( v ) = 0 i f f v = 0
如果拓扑向量空间的拓扑可以被范数导出，这个拓扑向量空间被称为赋范向量空间.
l 1 − n o r m : ∥ x ⃗ ∥ 1 = ∑ i = 1 n ∣ x i ∣ . l_1-norm:\quad \|\vec x\|_1=\sum_{i=1}^n |x_i|.
l 1 − n o r m : ∥ x ∥ 1 = i = 1 ∑ n ∣ x i ∣ .
l 2 − n o r m : ∥ x ⃗ ∥ 2 = ( ∑ i = 1 n ∣ x i ∣ 2 ) 1 2 = x T x . l_2-norm:\quad \|\vec x\|_2=(\sum_{i=1}^n |x_i|^2)^\frac{1}{2}=\sqrt{x^Tx}.
l 2 − n o r m : ∥ x ∥ 2 = ( i = 1 ∑ n ∣ x i ∣ 2 ) 2 1 = x T x .
l ∞ − n o r m : ∥ x ⃗ ∥ ∞ = max 1 ≤ i ≤ n ∣ x i ∣ . l_\infty-norm:\quad \|\vec x\|_\infty=\max_{1\leq i\leq n} |x_i|.
l ∞ − n o r m : ∥ x ∥ ∞ = 1 ≤ i ≤ n max ∣ x i ∣ .
l p − n o r m : ∥ x ⃗ ∥ p = ( ∑ i = 1 n ∣ x i ∣ p ) 1 p . l_p-norm:\quad \|\vec x\|_p=(\sum_{i=1}^n |x_i|^p)^\frac{1}{p}.
l p − n o r m : ∥ x ∥ p = ( i = 1 ∑ n ∣ x i ∣ p ) p 1 . matrix H is positive definite. H ≥ 0 \mathbf H\geq 0 H ≥ 0
H − n o r m : ∥ x ∥ H = x T ⋅ H ⋅ x . H-norm:\quad \|x\|_H=\sqrt {x^T\cdot H\cdot x}.
H − n o r m : ∥ x ∥ H = x T ⋅ H ⋅ x .
∥ A ∥ = max x ⃗ ≠ 0 ⃗ ∥ A x ⃗ ∥ ∥ x ⃗ ∥ \|\mathbf A\|=\max_{\vec x\neq\vec 0}\frac{\|\mathbf A\vec x\|}{\|\vec x\|}
∥ A ∥ = x = 0 max ∥ x ∥ ∥ A x ∥
1 − n o r m : ∥ A ∥ 1 = max j ∑ i = 1 m ∣ a i j ∣ . 1-norm:\quad \|\mathbf A\|_1=\max_{j}\sum_{i=1}^m|a_{ij}|.
1 − n o r m : ∥ A ∥ 1 = j max i = 1 ∑ m ∣ a i j ∣ .
∞ − n o r m : ∥ A ∥ ∞ = max i ∑ j = 1 m ∣ a i j ∣ . \infty-norm:\quad \|\mathbf A\|_\infty=\max_{i}\sum_{j=1}^m|a_{ij}|.
∞ − n o r m : ∥ A ∥ ∞ = i max j = 1 ∑ m ∣ a i j ∣ . Frobenius 范数定义为
F r o b e n i u s − n o r m : ∥ A ∥ F = ( ∑ i = 1 m ∑ j = 1 n a i j 2 ) 1 / 2 . Frobenius-norm:\quad \|\mathbf A\|_F= (\sum_{i=1}^m\sum_{j=1}^n a_{ij}^2)^{1/2}.
F r o b e n i u s − n o r m : ∥ A ∥ F = ( i = 1 ∑ m j = 1 ∑ n a i j 2 ) 1 / 2 . 已知实对称矩阵有完备的特征向量系，即有A V = V Λ AV=V\Lambda A V = V Λ V V V Λ \Lambda Λ
线性方程组左乘一个非奇异矩阵，方程组的解不变。
排列矩阵 Permutation Matrix
对角线调整 Postive Definite Diagonal Matrix
下三角线性方程组 L x ⃗ = b ⃗ \mathbf L\vec x=\vec b L x = b
上三角线性方程组 U x ⃗ = b ⃗ \mathbf U\vec x=\vec b U x = b
M k α = ( 1 ⋱ 1 − m k + 1 1 ⋮ ⋱ − m n 1 ) ( α 1 ⋮ α k α k + 1 ⋮ α n ) = ( α 1 ⋮ α k 0 ⋮ 0 ) M_k\alpha=
\left(
\begin{array}{ccc}
1 \\
& \ddots \\
& & 1\\
& & -m_{k+1} & 1\\
& & \vdots & & \ddots\\
& & -m_n & & & 1
\end{array}\right)
\left(
\begin{array}{ccc}
\alpha_1\\
\vdots\\
\alpha_k\\
\alpha_{k+1}\\
\vdots\\
\alpha_n\\
\end{array}\right)=
\left(
\begin{array}{ccc}
\alpha_1\\
\vdots\\
\alpha_k\\
0\\
\vdots\\
0\\
\end{array}\right)
M k α = ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎛ 1 ⋱ 1 − m k + 1 ⋮ − m n 1 ⋱ 1 ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎞ ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎛ α 1 ⋮ α k α k + 1 ⋮ α n ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎞ = ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎛ α 1 ⋮ α k 0 ⋮ 0 ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎞ M k M_k M k
M k = I − m k e k T = I − ( 0 ⋮ 0 0 m k + 1 ⋮ m n ) ( 0 , ⋯ , 0 , 1 k , 0 , ⋯ , 0 ) M_k=I-m_ke_k^T=I-
\left(\begin{array}{ccc}
0\\
\vdots\\
0\\
0\\
m_{k+1}\\
\vdots\\
m_n
\end{array}\right)
\left(\begin{array}{ccc}
0,\cdots,0,1_k,0,\cdots,0
\end{array}\right)
M k = I − m k e k T = I − ⎝ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎛ 0 ⋮ 0 0 m k + 1 ⋮ m n ⎠ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎞ ( 0 , ⋯ , 0 , 1 k , 0 , ⋯ , 0 ) 因为e k T m k = 0 e_k^Tm_k=0 e k T m k = 0 M K − 1 = I + m k e k T M_K^{-1}=I+m_ke_k^T M K − 1 = I + m k e k T e k T m j = 0 e_k^Tm_j=0 e k T m j = 0
M k M j = [ I − m k e k T ] [ I − m j e j T ] = I − m k e k T − m j e j T + m k ( e k T m j ) e j T = I − m k e k T − m j e j T M_kM_j=[I-m_ke_k^T][I-m_je_j^T]=I-m_ke_k^T-m_je_j^T+m_k(e_k^Tm_j)e_j^T=I-m_ke_k^T-m_je_j^T
M k M j = [ I − m k e k T ] [ I − m j e j T ] = I − m k e k T − m j e j T + m k ( e k T m j ) e j T = I − m k e k T − m j e j T **对应问题：** 求解$Ax=b$
M ⋅ A x = M n − 1 ⋯ M 2 M 1 A x = U x M\cdot Ax=M_{n-1}\cdots M_2M_1Ax=Ux
M ⋅ A x = M n − 1 ⋯ M 2 M 1 A x = U x
M ⋅ b = M n − 1 ⋯ M 2 M 1 b M\cdot b=M_{n-1}\cdots M_2M_1b
M ⋅ b = M n − 1 ⋯ M 2 M 1 b 高斯消去法即选主元（有限精度必须，避免误差放大），再将A转换成上/下三角矩阵后求解。
**对应问题：** 求解$Ax=b$, where A is an nxn matrix, symmetric and positive definite (A^T=A, A>0).
该问题可以等价于求函数f(x)的极值： f ( x ) = 1 2 x T A x − b T x f(x)=\frac{1}{2}x^TAx-b^Tx f ( x ) = 2 1 x T A x − b T x r ( x ) : = ∇ f ( x ) = A x − b r(x):=\nabla f(x)=Ax-b r ( x ) : = ∇ f ( x ) = A x − b
需要了解迭代、搜索方向等内容，此处不赘述。
如果非零向量集合{ p 0 , p 1 , ⋯ , p l } \{p_0,p_1,\cdots,p_l\} { p 0 , p 1 , ⋯ , p l } p i T A p j = 0 f o r a l l i ≠ j p_i^TAp_j=0\quad for\quad all\quad i\neq j p i T A p j = 0 f o r a l l i = j p i p_i p i
***Proof:***
对于$a_0p_0+a_1p_1+\cdots+a_lp_l=0$，左乘$p_0^T$得$p_0^Ta_0p_0+\underbrace{\cdots}_{=0}=0$。又因为$p_0^Ta_0p_0=a_0(p_0^Tp_0)\quad and\quad (p_0^Tp_0)\neq 0$ 所以 $a_0=0$。同理可得$\forall a_i=0$，故得证。
用共轭向量p k p_k p k x k + 1 = x k + α k p k x_{k+1}=x_k+\alpha_kp_k x k + 1 = x k + α k p k α k ( : = − r k T p k p k T A p k ) \alpha_k(:=-\frac{r_k^Tp_k}{p_k^TAp_k}) α k ( : = − p k T A p k r k T p k )
共轭梯度算法可在最多n次迭代内收敛到问题的目标解。
For any $x_0 \in R^n$, the sequence {$x_k$} generated by the conjugate direction algorithm converges to the solution $x^*$ of the linear system in **at most n steps**.
***Proof:***
当迭代n次时，因为{$p_i$}线性无关，可以作为n维空间的基，张成空间内的任意向量，包括目标解$x^*$。故若选取合适的$\alpha_i$满足$x_k-x^*=(x_0+\sum_{i=1}^{k-1}\alpha_ip_i)-(\sum_{i=0}^{n-1}\delta_i^*p_i)=0$，最多n次迭代即可得到目标解。
已知存在表达式x ∗ − x 0 = ∑ i = 0 n − 1 δ i p i x^*-x_0=\sum_{i=0}^{n-1}\delta_ip_i x ∗ − x 0 = ∑ i = 0 n − 1 δ i p i p k T A p_k^TA p k T A
p k T A ( x ∗ − x 0 ) = p k T A ∑ i = 0 n − 1 δ i p i = p k T A δ k p k p k T A ( x ∗ − x 0 ) = p k T A ( x ∗ − x k + x k − x 0 ) = p k T A ( x ∗ − x k ) + 0 = p k T ( b − A x k ) = − p k T r k L e t α k = δ k = − r k T p k p k T A p k p_k^TA(x^*-x_0)=p_k^TA\sum_{i=0}^{n-1}\delta_ip_i=p_k^TA\delta_kp_k\\
p_k^TA(x^*-x_0)=p_k^TA(x^*-x_k+x_k-x_0)=p_k^TA(x^*-x_k)+0=p_k^T(b-Ax_k)=-p_k^Tr_k\\
Let\quad \alpha_k=\delta_k=\frac{-r_k^Tp_k}{p_k^TAp_k}
p k T A ( x ∗ − x 0 ) = p k T A i = 0 ∑ n − 1 δ i p i = p k T A δ k p k p k T A ( x ∗ − x 0 ) = p k T A ( x ∗ − x k + x k − x 0 ) = p k T A ( x ∗ − x k ) + 0 = p k T ( b − A x k ) = − p k T r k L e t α k = δ k = p k T A p k − r k T p k x 0 ∈ R n x_0\in \mathbb{R}^n x 0 ∈ R n { x i } i = 1 k \{x_i\}_{i=1}^k { x i } i = 1 k
part (1)
r k T p i = 0 , f o r i = 0 , 1 , … , k − 1 r_k^Tp_i=0,\quad for\quad i=0, 1,\dots,k-1
r k T p i = 0 , f o r i = 0 , 1 , … , k − 1 part (2)
对于{ x ∣ x = x 0 + s p a n ( p 0 , p 1 , … , p k − 1 ) } \{x|x=x_0+span(p_0,p_1,\dots,p_{k-1})\} { x ∣ x = x 0 + s p a n ( p 0 , p 1 , … , p k − 1 ) } x k x_k x k f ( x ) = 1 2 x T A x − b T x f(x)=\frac{1}{2}x^TAx-b^Tx f ( x ) = 2 1 x T A x − b T x
***Proof:***
(1)显然的，有以下几个式子成立
r 1 T p 0 = ( r 0 + α 0 A p 0 ) T p 0 = r 0 T p 0 − r 0 T p 0 = 0 r_1^Tp_0=(r_0+\alpha_0 Ap_0)^Tp_0=r_0^Tp_0-r_0^Tp_0=0
r 1 T p 0 = ( r 0 + α 0 A p 0 ) T p 0 = r 0 T p 0 − r 0 T p 0 = 0
r k T p k − 1 = ( r k − 1 + α k − 1 A p k − 1 ) T p k − 1 = 0 r_k^Tp_{k-1}=(r_{k-1}+\alpha_{k-1} Ap_{k-1})^Tp_{k-1}=0
r k T p k − 1 = ( r k − 1 + α k − 1 A p k − 1 ) T p k − 1 = 0
r k T p i = ( r k − 1 + α k − 1 A p k − 1 ) T p i = r k − 1 T p i − α k − 1 p k − 1 T A p i = r k − 1 T p i + 0 f o r i = 0 , 1 , … , k − 2 r_k^Tp_i=(r_{k-1}+\alpha_{k-1} Ap_{k-1})^Tp_i\\
=r_{k-1}^Tp_i-\alpha_{k-1} p_{k-1}^TAp_i\\
=r_{k-1}^Tp_i+0\\
for\quad i=0,1,\dots,k-2 r k T p i = ( r k − 1 + α k − 1 A p k − 1 ) T p i = r k − 1 T p i − α k − 1 p k − 1 T A p i = r k − 1 T p i + 0 f o r i = 0 , 1 , … , k − 2 最后一个式子说明r k T p i r_k^Tp_i r k T p i r i + 1 T p i = 0 r_{i+1}^Tp_i=0 r i + 1 T p i = 0
(2)令h ( σ ) = f ( x 0 + σ 0 p 0 + ⋯ + σ k − 1 p k − 1 ) h(\sigma)=f(x_0+\sigma_0p_0+\dots+\sigma_{k-1}p_{k-1}) h ( σ ) = f ( x 0 + σ 0 p 0 + ⋯ + σ k − 1 p k − 1 ) σ ∗ \sigma^* σ ∗
∂ h ( σ ∗ ) ∂ σ i = 0 = ∂ h ( σ ∗ ) ∂ ( σ i p i ) ∂ ( σ i p i ) ∂ σ i = ∂ f ( x 0 + σ 0 p 0 + ⋯ + σ k − 1 p k − 1 ) ∂ ( σ i p i ) p i = ∇ f ( x ∗ ) T p i = ( A x ∗ − b ) T p i = r ( x ∗ ) T p i f o r i = 0 , 1 , … , k − 1 \frac{\partial h(\sigma^*)}{\partial\sigma_i}=0=\frac{\partial h(\sigma^*)}{\partial(\sigma_ip_i)}\frac{\partial(\sigma_ip_i)}{\partial\sigma_i}\\
=\frac{\partial f(x_0+\sigma_0p_0+\dots+\sigma_{k-1}p_{k-1})}{\partial(\sigma_ip_i)}p_i\\
=\nabla f(x^*)^Tp_i=(Ax^*-b)^Tp_i\\
=r(x^*)^Tp_i\\
for\quad i=0,1,\dots,k-1 ∂ σ i ∂ h ( σ ∗ ) = 0 = ∂ ( σ i p i ) ∂ h ( σ ∗ ) ∂ σ i ∂ ( σ i p i ) = ∂ ( σ i p i ) ∂ f ( x 0 + σ 0 p 0 + ⋯ + σ k − 1 p k − 1 ) p i = ∇ f ( x ∗ ) T p i = ( A x ∗ − b ) T p i = r ( x ∗ ) T p i f o r i = 0 , 1 , … , k − 1 同时有r k = r ( x k ) = A x k − b r_k=r(x_k)=Ax_k-b r k = r ( x k ) = A x k − b r k T p i = 0 r_k^Tp_i=0 r k T p i = 0 r ( x k ) T p i = r ( x ∗ ) T p i = 0 r(x_k)^Tp_i=r(x^*)^Tp_i=0 r ( x k ) T p i = r ( x ∗ ) T p i = 0
故x k x_k x k x ∗ x^* x ∗
由于{ p i } i = 0 k \{p_i\}_{i=0}^k { p i } i = 0 k p k p_k p k p k − 1 p_{k-1} p k − 1 { p i } i = 0 k − 2 \{p_i\}_{i=0}^{k-2} { p i } i = 0 k − 2
取p k = − r k + β k p k − 1 p_k=-r_k+\beta_k p_{k-1} p k = − r k + β k p k − 1 p k − 1 T A p_{k-1}^TA p k − 1 T A
β k = p k − 1 T A r k p k − 1 T A p k − 1 \beta_k=\frac{p_{k-1}^TAr_k}{p_{k-1}^TAp_{k-1}}
β k = p k − 1 T A p k − 1 p k − 1 T A r k 因此有以下性质：
part (1)
r k T r i = 0 f o r i = 0 , … , k − 1 r_k^Tr_i=0\quad for\quad i=0,\dots,k-1
r k T r i = 0 f o r i = 0 , … , k − 1 part (2)
s p a n ( p 0 , … , p k ) = s p a n ( r 0 , … , r k ) = s p a n ( r 0 , A r 0 , … , A k r 0 ) span(p_0,\dots,p_k)=span(r_0,\dots,r_k)=span(r_0,Ar_0,\dots,A^kr_0)
s p a n ( p 0 , … , p k ) = s p a n ( r 0 , … , r k ) = s p a n ( r 0 , A r 0 , … , A k r 0 ) ***Proof:***
首先$r_k=\beta_kp_{k-1}-p_k$ 所以仅有i=k-1时需要证明其不为0。
$$r_k^Tr_{k-1}=\beta_kp_{k-1}^Tr_{k-1}=\frac{p_{k-1}^TA(r_kp_{k-1}^T)r_{k-1}}{p_{k-1}^TAp_{k-1}}=0$$
(1)得证。
用数学归纳法易证(2)。假设r k ∈ s p a n ( r 0 , A r 0 , … , A k r 0 ) r_k\in span(r_0,Ar_0,\dots,A^kr_0) r k ∈ s p a n ( r 0 , A r 0 , … , A k r 0 ) r k + 1 r_{k+1} r k + 1 s p a n ( r 0 , … , r k ) ∈ s p a n ( r 0 , A r 0 , … , A k r 0 ) span(r_0,\dots,r_k)\in span(r_0,Ar_0,\dots,A^kr_0) s p a n ( r 0 , … , r k ) ∈ s p a n ( r 0 , A r 0 , … , A k r 0 )
CG-Preliminary Version
Set r 0 : = A x 0 − b , p 0 : = − r 0 , k : = 0 r_0:=Ax_0-b,\quad p_0:=-r_0,\quad k:=0 r 0 : = A x 0 − b , p 0 : = − r 0 , k : = 0
while r k ≥ t o l e r a n c e r_k\geq tolerance r k ≥ t o l e r a n c e
α k : = − r k T p k p k T A p k \alpha_k:=-\frac{r_k^Tp_k}{p_k^TAp_k} α k : = − p k T A p k r k T p k x k + 1 : = x k + α k p k x_{k+1}:=x_k+\alpha_kp_k x k + 1 : = x k + α k p k r k + 1 : = A x k + 1 − b r_{k+1}:=Ax_{k+1}-b r k + 1 : = A x k + 1 − b β k + 1 : = r k + 1 T A p k p k T A p k \beta_{k+1}:=\frac{r_{k+1}^TAp_k}{p_k^TAp_k} β k + 1 : = p k T A p k r k + 1 T A p k p k + 1 : = − r k + 1 + β k + 1 p k p_{k+1}:=-r_{k+1}+\beta_{k+1}p_k p k + 1 : = − r k + 1 + β k + 1 p k k + + k++ k + +
end CG-Pratical Version
Set r 0 : = A x 0 − b , p 0 : = − r 0 , k : = 0 r_0:=Ax_0-b,\quad p_0:=-r_0,\quad k:=0 r 0 : = A x 0 − b , p 0 : = − r 0 , k : = 0
while r k ≥ t o l e r a n c e r_k\geq tolerance r k ≥ t o l e r a n c e
α k : = − r k T r k p k T A p k ( = − r k T p k p k T A p k ) \alpha_k:=-\frac{r_k^Tr_k}{p_k^TAp_k}(=-\frac{r_k^Tp_k}{p_k^TAp_k}) α k : = − p k T A p k r k T r k ( = − p k T A p k r k T p k ) x k + 1 : = x k + α k p k x_{k+1}:=x_k+\alpha_kp_k x k + 1 : = x k + α k p k r k + 1 : = r k + α k A p k r_{k+1}:=r_k+\alpha_kAp_k r k + 1 : = r k + α k A p k β k + 1 : = r k + 1 T r k + 1 r k T r k \beta_{k+1}:=\frac{r_{k+1}^Tr_{k+1}}{r_k^Tr_k} β k + 1 : = r k T r k r k + 1 T r k + 1 p k + 1 : = − r k + 1 + β k + 1 p k p_{k+1}:=-r_{k+1}+\beta_{k+1}p_k p k + 1 : = − r k + 1 + β k + 1 p k k + + k++ k + +
end 易证两者是等价的。
**对应问题：** 求解$Ax=b.$ A为非奇异矩阵。
当A为低阶稠密矩阵时，高斯消去法可以有效求解。但对于大型稀疏矩阵（A阶数很大，且零元素较多，比如n>10^4），这时利用迭代法求解时合适的。在计算机内存和运算两方面，迭代法都可以利用大量零元素的特点。
待更新
待更新
待更新
泰勒展开
一阶必要条件：∇ f ( x ∗ ) = 0 \nabla f(x^*)=0 ∇ f ( x ∗ ) = 0
***Proof:*** 用反证法，利用泰勒公式构造式子，证明与局部最小矛盾。
二阶必要条件：∇ 2 f ( x ) \nabla^2f(x) ∇ 2 f ( x ) ∇ f ( x ∗ ) = 0 a n d ∇ 2 f ( x ∗ ) \nabla f(x^*)=0\quad and\quad \nabla^2f(x^*) ∇ f ( x ∗ ) = 0 a n d ∇ 2 f ( x ∗ )
***Proof:*** 同理，用反证法，利用泰勒公式构造式子，证明与局部最小矛盾。
二阶充分条件：∇ 2 f ( x ) \nabla^2f(x) ∇ 2 f ( x ) ∇ f ( x ∗ ) = 0 a n d ∇ 2 f ( x ∗ ) \nabla f(x^*)=0\quad and\quad \nabla^2f(x^*) ∇ f ( x ∗ ) = 0 a n d ∇ 2 f ( x ∗ )
***Proof:*** 取半径r>0构造球形域，此时Hessian矩阵$\nabla^2f(x)$对球形域内任意x而言都是正定的。对于任意范数$\|p\|$小于r的向量，有
f(x^*+p)=f(x^*)+p^T\nabla f(x^*)+\frac{1}{2}p^T\nabla f(x^*+tp)p\\
=f(x^*)+\frac{1}{2}p^T\nabla f(x^*+tp)p$$。
而因为Hessian矩阵正定，有$p^T\nabla f(x^*+tp)p>0$，故局部最小得证。
5. (A). 如果函数convex，那么任意局部最小x^star就是全局最小；(B). 如果函数可微，那么任意驻点x^star是全局最小
***Proof:*** 利用Convexity property 和 Jensen's inequality进行构造式子可证(A)；用反证法，假设$x^*$是驻点但不是全局最小。首先$\lim_{\lambda\to 0}\frac{f(x^*+\lambda(z-x^*))-f(x^*)}{\lambda}=\nabla f(x^*)^T(z-x^*)$。
利用Jensen's inequality得
$$\lim_{\lambda\to 0}\frac{f(x^*+\lambda(z-x^*))-f(x^*)}{\lambda}\leq\lim_{\lambda\to 0}\frac{\lambda f(z)+(1-\lambda)f(x^*)-f(x^*)}{\lambda}\\=f(z)-f(x^*)<0$$，得到$\nabla f(x^*)^T<0$即$x^*$不是驻点，与假设矛盾，故(B)得证。
## 下降方法的收敛性
### 关于收敛
线性收敛 Linear Convergence
$$\|x_{k+1}-x^*\|\leq c\|x_k-x^*\|\quad 0
二次收敛 Quadratic Convergence
∥ x k + 1 − x ∗ ∥ ≤ c ∥ x k − x ∗ ∥ 2 0 < c < 1 \|x_{k+1}-x^*\|\leq c\|x_k-x^*\|^2\quad 0<c<1
∥ x k + 1 − x ∗ ∥ ≤ c ∥ x k − x ∗ ∥ 2 0 < c < 1 超线性收敛 Superlinear Convergence
lim k → 0 ∥ x k + 1 − x ∗ ∥ ∥ x k − x ∗ ∥ = 0 \lim_{k\to 0}\frac{\|x_{k+1}-x^*\|}{\|x_k-x^*\|}=0
k → 0 lim ∥ x k − x ∗ ∥ ∥ x k + 1 − x ∗ ∥ = 0
Method
Complexity
Convergence
Steepest Descent
O(n)
Linear
Conjugate Gradient
O(n)
Linearfail
Newton’s
O(n^3)
Quadratic
Damped Newton’s
O(n^3)
Superlinear
Quasi-Newton (DFP)
O(n^2)
Superlinear
Quasi-Newton (BFGS)
O(n^2)
Superlinear
此处内容可参考 ***Convex Optimization***, *Stephen Boyd (page475)*
最速下降法比较简单，只需搜索方向设置为最速下降方向（假设读者知道优化中的搜索相关内容）。对于f ( x + v ) ≃ f ( x ) + ∇ f ( x ) T v f(x+v)\simeq f(x)+\nabla f(x)^Tv f ( x + v ) ≃ f ( x ) + ∇ f ( x ) T v
Δ x n s d = a r g m i n { ∇ f ( x ) T v ∣ ∥ v ∥ = 1 } \Delta x_{nsd}=argmin\{\nabla f(x)^Tv\Big|\|v\|=1\}
Δ x n s d = a r g m i n { ∇ f ( x ) T v ∣ ∣ ∣ ∥ v ∥ = 1 } 当v的方向为梯度方向时，即是梯度下降法。而相比而言，最速下降可以选取非Euclidean norm，或者多个最优解中的其中一个。
Reference: 最速下降法与梯度下降法区别
**对应问题：** 求解$g(x)=0.$ 以及可以转化成该形式的问题。
此处内容同计算统计 ***Computational Statistics***, *Guo-Liang TIAN, 2018 (page33-37)*
设x t x_t x t x ∗ x^* x ∗
g ( x t ) − 0 = g ( x t ) − g ( x ∗ ) = g ′ ( x α ) ( x t − x ∗ ) x α ∈ [ m i n { x t , x ∗ } , m a x { x t , x ∗ } ] g(x_t)-0=g(x_t)-g(x^*)=g'(x_\alpha)(x_t-x^*)\\
x_\alpha \in \big[min\{x_t, x^*\}, max\{x_t, x^*\}\big]
g ( x t ) − 0 = g ( x t ) − g ( x ∗ ) = g ′ ( x α ) ( x t − x ∗ ) x α ∈ [ m i n { x t , x ∗ } , m a x { x t , x ∗ } ] 用x t x_t x t x α x_\alpha x α x t + 1 x_{t+1} x t + 1 x ∗ x^* x ∗
x t + 1 = x t − g ( x t ) g ′ ( x t ) x_{t+1}=x_t-\frac{g(x_t)}{g'(x_t)}
x t + 1 = x t − g ′ ( x t ) g ( x t )
L e t F ( x ) = x − g ( x ) g ′ ( x ) Let\quad F(x)=x-\frac{g(x)}{g'(x)}
L e t F ( x ) = x − g ′ ( x ) g ( x )
F ′ ( x ∗ ) = F ′ ( x ) ∣ x = x ∗ = 1 − g ′ ( x ) 2 − g ( x ) g ′ ′ ( x ) g ′ ( x ) 2 ∣ x = x ∗ = g ( x ) g ′ ′ ( x ) g ′ ( x ) 2 ∣ x = x ∗ = g ( x ∗ ) g ′ ′ ( x ∗ ) g ′ ( x ∗ ) 2 = 0 ∵ g ( x ∗ ) = 0 F'(x^*)=F'(x)|_{x=x^*}=1-\frac{g'(x)^2-g(x)g''(x)}{g'(x)^2}\bigg|_{x=x^*}\\
=\frac{g(x)g''(x)}{g'(x)^2}\bigg|_{x=x^*}=\frac{g(x^*)g''(x^*)}{g'(x^*)^2}\\
=0\qquad \because g(x^*)=0
F ′ ( x ∗ ) = F ′ ( x ) ∣ x = x ∗ = 1 − g ′ ( x ) 2 g ′ ( x ) 2 − g ( x ) g ′ ′ ( x ) ∣ ∣ ∣ ∣ x = x ∗ = g ′ ( x ) 2 g ( x ) g ′ ′ ( x ) ∣ ∣ ∣ ∣ x = x ∗ = g ′ ( x ∗ ) 2 g ( x ∗ ) g ′ ′ ( x ∗ ) = 0 ∵ g ( x ∗ ) = 0 收敛率 c = lim t → ∞ ∥ x t − x ∗ ∥ ∥ x t − x ∗ ∥ 2 = 1 2 ∣ f ′ ′ ( x ^ ) ∣ c=\lim_{t\to\infty} \frac{\|x_t-x^*\|}{\|x_t-x^*\|^2}=\frac{1}{2}|f''(\hat x)| c = lim t → ∞ ∥ x t − x ∗ ∥ 2 ∥ x t − x ∗ ∥ = 2 1 ∣ f ′ ′ ( x ^ ) ∣
该部分内容可以参考 ***Numerical Optimization***, *1999, Nocedal J., Wright S.J. (page134-162)*
包括
Inexact Newton Steps
Line Search Newton Methods
Hessian Modifications
Trust-Region Newton Methods
用对称矩阵B k = ∇ 2 f ( x k ) B_k=\nabla^2 f(x_k) B k = ∇ 2 f ( x k )
p k = a r g m i n { m k ( p ) = f ( x k ) + f ( x k ) T p + 1 2 p T B k p } s . t . ∥ p ∥ ≤ Δ k p_k=argmin\{m_k(p)=f(x_k)+f(x_k)^Tp+\frac{1}{2}p^TB_kp\}\\
s.t.\quad \|p\|\leq\Delta_k
p k = a r g m i n { m k ( p ) = f ( x k ) + f ( x k ) T p + 2 1 p T B k p } s . t . ∥ p ∥ ≤ Δ k 根据比例ρ k = f ( x k ) − f ( x k + p k ) m k ( 0 ) − m k ( p k ) \rho_k=\frac{f(x_k)-f(x_k+p_k)}{m_k(0)-m_k(p_k)} ρ k = m k ( 0 ) − m k ( p k ) f ( x k ) − f ( x k + p k ) Δ k \Delta_k Δ k x k + 1 = x k + p k . x_{k+1}=x_k+p_k. x k + 1 = x k + p k .
其中有个the dogleg method：狗腿法，名字比较有趣 😃
Reference: 数值优化 Ch.4 信赖域方法
**对应问题：** 求解$g(x)=0.$ 以及可以转化成该形式的问题。
牛顿法虽然收敛速度快，但每次迭代都要计算Hessian矩阵，极大地复杂了计算度。所以类似Fisher Scoring算法对Newton-Raphson算法的改进，拟牛顿法构造出可以近似Hessian矩阵的正定对称阵，在“拟牛顿”的条件下优化目标函数。
下文B表示对Hessian矩阵的近似，D表示Hessian矩阵的逆矩阵的近似。即 B ≃ H , D ≃ H − 1 B\simeq H, D\simeq H^{-1} B ≃ H , D ≃ H − 1
k+1次迭代后进行泰勒展开有：
∇ f ( x ) ≃ ∇ f ( x k + 1 ) + H k + 1 ⋅ ( x − x k + 1 ) \nabla f(x)\simeq\nabla f(x_{k+1})+H_{k+1}\cdot (x-x_{k+1})
∇ f ( x ) ≃ ∇ f ( x k + 1 ) + H k + 1 ⋅ ( x − x k + 1 ) 取x = x k x=x_k x = x k s k = x k + 1 − x k , y k = ∇ f ( x k + 1 ) − ∇ f ( x k ) s_k=x_{k+1}-x_k, y_k=\nabla f(x_{k+1})-\nabla f(x_k) s k = x k + 1 − x k , y k = ∇ f ( x k + 1 ) − ∇ f ( x k )
y k ≃ H k + 1 ⋅ s k y_k\simeq H_{k+1}\cdot s_k
y k ≃ H k + 1 ⋅ s k 或
s k ≃ H k + 1 − 1 ⋅ y k s_k\simeq H_{k+1}^{-1}\cdot y_k
s k ≃ H k + 1 − 1 ⋅ y k 这就是拟牛顿条件，它对迭代过程中的Hessian矩阵作约束。故近H似矩阵B和逆矩阵近似矩阵D满足：
y k ≃ B k + 1 ⋅ s k y_k\simeq B_{k+1}\cdot s_k
y k ≃ B k + 1 ⋅ s k 或
s k ≃ D k + 1 − 1 ⋅ y k s_k\simeq D_{k+1}^{-1}\cdot y_k
s k ≃ D k + 1 − 1 ⋅ y k DFP核心在于：用D k D_k D k H k − 1 H_{k}^{-1} H k − 1
D k + 1 = D k + Δ D k , k = 0 , 1 , 2 , ⋯ D_{k+1}=D_k+\Delta D_k,\quad k=0,1,2,\cdots
D k + 1 = D k + Δ D k , k = 0 , 1 , 2 , ⋯ 为保持D的对称性，假设有：
Δ D k = α u u T + β v v T \Delta D_k=\alpha uu^T+\beta vv^T
Δ D k = α u u T + β v v T 将其代入有
s k = D k y k + α u u T y k + β v v T y k = D k y k + u ( α u T y k ) + v ( β v T y k ) = D k y k + ( α u T y k ) u + ( β v T y k ) v s_k=D_ky_k+\alpha uu^Ty_k+\beta vv^Ty_k\\
=D_ky_k+u(\alpha u^Ty_k)+v(\beta v^Ty_k)\\
=D_ky_k+(\alpha u^Ty_k)u+(\beta v^Ty_k)v
s k = D k y k + α u u T y k + β v v T y k = D k y k + u ( α u T y k ) + v ( β v T y k ) = D k y k + ( α u T y k ) u + ( β v T y k ) v 不妨令α u T y k = 1 , β v T y k = − 1 \alpha u^Ty_k=1,\quad\beta v^Ty_k=-1 α u T y k = 1 , β v T y k = − 1
u − v = s k − D k y k u-v=s_k-D_ky_k
u − v = s k − D k y k 不妨令u = s k , v = D k y k u=s_k,\quad v=D_ky_k u = s k , v = D k y k
α = 1 s k T y k , β = − 1 ( D k y k ) T y k = − 1 y k T D k T y k = − 1 y k T D k y k \alpha=\frac{1}{s_k^Ty_k},\quad\beta=-\frac{1}{(D_ky_k)^Ty_k}=\frac{-1}{y_k^TD_k^Ty_k}=\frac{-1}{y_k^TD_ky_k}
α = s k T y k 1 , β = − ( D k y k ) T y k 1 = y k T D k T y k − 1 = y k T D k y k − 1 这样就构造出了一个校正矩阵Δ D k = s k s k T s k T y k − D k y k y k T D k y k T D k y k \Delta D_k=\frac{s_ks_k^T}{s_k^Ty_k}-\frac{D_ky_ky_k^TD_k}{y_k^T D_k y_k} Δ D k = s k T y k s k s k T − y k T D k y k D k y k y k T D k
DFP
We have: D k ≃ H k − 1 , g k = ∇ f ( x k ) D_k\simeq H_k^{-1}, g_k=\nabla f(x_k) D k ≃ H k − 1 , g k = ∇ f ( x k )
Set x_0, tolerance and D 0 : = I , k : = 0 D_0:=I,\quad k:=0 D 0 : = I , k : = 0
while g k ≥ t o l e r a n c e g_k\geq tolerance g k ≥ t o l e r a n c e
d k : = − D k ⋅ g k d_k:=-D_k\cdot g_k d k : = − D k ⋅ g k λ k : = a r g m i n { f ( x k + λ d k ) } \lambda_k:=argmin\{f(x_k+\lambda d_k)\} λ k : = a r g m i n { f ( x k + λ d k ) } s k : = λ k d k s_k:=\lambda_kd_k s k : = λ k d k x k + 1 : = x k + s k x_{k+1}:=x_k+s_k x k + 1 : = x k + s k y k : = g k + 1 − g k y_k:=g_{k+1}-g_k y k : = g k + 1 − g k D k + 1 : = D k + s k s k T s k T y k − D k y k y k T D k y k T D k y k D_{k+1}:=D_k+\frac{s_ks_k^T}{s_k^Ty_k}-\frac{D_ky_ky_k^TD_k}{y_k^T D_k y_k} D k + 1 : = D k + s k T y k s k s k T − y k T D k y k D k y k y k T D k k + + k++ k + +
end BFGS核心在于直接逼近近似H k H_{k} H k
B k + 1 = B k + Δ B k , k = 0 , 1 , 2 , ⋯ B_{k+1}=B_k+\Delta B_k,\quad k=0,1,2,\cdots
B k + 1 = B k + Δ B k , k = 0 , 1 , 2 , ⋯ 同样，为保持B的对称性，假设有：
Δ B k = α u u T + β v v T \Delta B_k=\alpha uu^T+\beta vv^T
Δ B k = α u u T + β v v T 将其代入有
y k = B k s k + ( α u T s k ) u + ( β v T s k ) v y_k=B_ks_k+(\alpha u^Ts_k)u+(\beta v^Ts_k)v
y k = B k s k + ( α u T s k ) u + ( β v T s k ) v 不妨令α u T s k = 1 , β v T s k = − 1 \alpha u^Ts_k=1,\quad\beta v^Ts_k=-1 α u T s k = 1 , β v T s k = − 1 u = y k , v = B k s k u=y_k,\quad v=B_ks_k u = y k , v = B k s k
α = 1 y k T s k , β = − 1 s k T B k s k \alpha=\frac{1}{y_k^Ts_k},\quad\beta=\frac{-1}{s_k^TB_ks_k}
α = y k T s k 1 , β = s k T B k s k − 1 这样就构造出了一个校正矩阵Δ B k = y k y k T y k T s k − B k s k s k T B k s k T B k s k \Delta B_k=\frac{y_ky_k^T}{y_k^Ts_k}-\frac{B_ks_ks_k^TB_k}{s_k^T B_k s_k} Δ B k = y k T s k y k y k T − s k T B k s k B k s k s k T B k
DFP
We have: B k ≃ H k , g k = ∇ f ( x k ) B_k\simeq H_k, g_k=\nabla f(x_k) B k ≃ H k , g k = ∇ f ( x k )
Set x_0, tolerance and B 0 : = I , k : = 0 B_0:=I,\quad k:=0 B 0 : = I , k : = 0
while g k ≥ t o l e r a n c e g_k\geq tolerance g k ≥ t o l e r a n c e
d k : = − B k − 1 ⋅ g k d_k:=-B_k^{-1}\cdot g_k d k : = − B k − 1 ⋅ g k λ k : = a r g m i n { f ( x k + λ d k ) } \lambda_k:=argmin\{f(x_k+\lambda d_k)\} λ k : = a r g m i n { f ( x k + λ d k ) } s k : = λ k d k s_k:=\lambda_kd_k s k : = λ k d k x k + 1 : = x k + s k x_{k+1}:=x_k+s_k x k + 1 : = x k + s k y k : = g k + 1 − g k y_k:=g_{k+1}-g_k y k : = g k + 1 − g k B k + 1 : = B k + y k y k T y k T s k − B k s k s k T B k s k T B k s k B_{k+1}:=B_k+\frac{y_ky_k^T}{y_k^Ts_k}-\frac{B_ks_ks_k^TB_k}{s_k^T B_k s_k} B k + 1 : = B k + y k T s k y k y k T − s k T B k s k B k s k s k T B k k + + k++ k + +
end 待更
