凸二次优化问题
Theory. 设
等价于求解线性方程组
Proof. 二次型二阶可导，极小值点处梯度为零，现对优化的目标函数求梯度。二次型本质上具有：
计算梯度的分量表达式：
合在一起写成矩阵形式：
显然，凸二次优化问题的极值条件等价于该线性方程组。凸二次优化问题在建模中十分常见，这说明讨论线性方程组的求解方法具有普遍的实用价值。然而对于规模较大的问题，使用线性代数中的克莱姆法则暴力展开将导致时间开销巨大，而高斯消元法算法流程又较为复杂。本文将介绍一种常见的数值分析方法，求得线性方程组的数值近似解。
最速梯度下降法
称优化目标函数的梯度为残量（Residue），即是当前解对于线性方程组的不满足量：
由于函数是凸函数，极值点一定存在。当前解处函数的梯度值表示了函数值上升最快的方向（梯度方向上方向导数最大）。那么沿着相反的方向每迭代一步就会更加靠近最优的极小值解。于是，我们定义迭代关系：
其中
在最速下降法中，将
同样的，驻点处梯度为零。根据链式求导法则：
最后一步由于
算例设计与实验
参考南科大的数值分析作业题EH计算设计SUSTech的算例一（由于版权关系无传送门）
生成实对称正定矩阵 采用文献[1]中的类似方法生成矩阵
其中
Householder矩阵是对称正交矩阵，这时，对角矩阵的特征值就是
推导条件数计算公式 矩阵
那么其特征值介于
于是可以反解处条件数的计算公式：
为了能够重复试验数据，使用线性同余法计算伪随机数进行向量的初始化（也可以使用相同的随机种子，调用编程语言内部的随机值库，这里尊重南科大的原题要求）。令xopt表示最优解
v = zeros(n, 1); xopt = zeros(n, 1); x0 = zeros(n, 1);
t = 0; for j = 1: n t = mod(t * 31416 + 13846, 46261); v(j) = t * (2 / 46261) - 1; end;
t = 0; for j = 1: n t = mod(t * 42108 + 13846, 46273); xopt(j) = t * (5 / 46273) + 5; end;
for j = 1: n t = mod(t * 42108 + 13846, 46273); x0(j) = t * (5 / 46273) - 10; end;这里设定最优解的分量在
这里设定算法的停机标准为当前梯度变为接近于零的极小量：
Python程序代码
import numpy as np
def initialize(cond, numb):
gamma = (np.cos(np.pi / (numb + 1)) + cond * np.cos(numb * np.pi / (numb + 1)) - (cond - 1)) / (cond - 1)
diags = np.array(list(map(lambda i: np.cos(i * np.pi / (numb + 1)) + 1 + gamma, np.arange(1, numb + 1))))
SIGMA = np.diag(diags)
v = np.zeros((numb, 1))
xopt = np.zeros((numb, 1))
x0 = np.zeros((numb, 1))
t = 0
for j in range(numb):
t = (t * 31416 + 13846) % 46261
v[j] = t * (2 / 46261) - 1
t = 0
for j in range(numb):
t = (t * 42108 + 13846) % 46273
xopt[j] = t * (5 / 46273) + 5
for j in range(numb):
t = (t * 42108 + 13846) % 46273
x0[j] = t * (5 / 46273) - 10
V = np.identity(numb) - 2 * v @ v.T / np.linalg.norm(v) ** 2
A = V @ SIGMA @ V.T
b = A @ xopt
return A, b, xopt, x0
def gradientDescent(A, b, x0, epsilon):
x_now = x0; epoch = 0
gnorm = ginit = np.linalg.norm(A @ x0 - b)
while gnorm / ginit > epsilon:
g_now = A @ x_now - b
gnorm = np.linalg.norm(g_now)
xnext = x_now - (gnorm ** 2) / (g_now.T @ A @ g_now) * g_now # * LR_linearDecay(epoch)
x_now = xnext
epoch += 1
return epoch, x_now
def LR_linearDecay(epoch):
return - epoch / 5e4 + 0.9
def LR_expertDecay(epoch):
return np.exp(- np.log(0.9) * (epoch / 4e3 - 1))
if __name__ == '__main__':
print("{:>8}".format("EPOCH/ERR"), end=" ")
for n in range(1, 6):
print("n={:>9}".format(100*n), end=" ")
print()
for c in range(3, 7):
print(f"COND=1e+{c}", end=" ")
for n in range(1, 6):
A, b, xopt, x0 = initialize(cond=10**c, numb=100*n)
epoch, xk = gradientDescent(A, b, x0, epsilon=1e-6)
error = np.linalg.norm(xk - xopt) / np.linalg.norm(x0 - xopt)
print("{:>5}/{:.4f}".format(epoch, error), end=" ")
print()MATLAB程序代码
Linear_init.m
function [A, b, xopt, x0] = linear_init(cond, numb)
gamma = (cos(pi / (numb + 1)) ...
+ cond * cos(numb * pi / (numb + 1)) - (cond - 1)) / (cond - 1);
diags = zeros(numb, 1);
for i = 1: numb;
diags(i) = cos(i * pi / (numb + 1)) + 1 + gamma;
end;
SIGMA = diag(diags);
v = zeros(numb, 1); xopt = zeros(numb, 1); x0 = zeros(numb, 1);
t1 = 0; t2 = 0;
for j = 1: numb
t1 = mod(t1 * 31416 + 13846, 46261);
v(j) = t1 * (2 / 46261) - 1;
end;
for j = 1: numb
t2 = mod(t2 * 42108 + 13846, 46273);
xopt(j) = t2 * (5 / 46273) + 5;
end;
for j = 1: numb
t2 = mod(t2 * 42108 + 13846, 46273);
x0(j) = t2 * (5 / 46273) - 10;
end;
V = eye(numb) - 2 * (v * v') / norm(v, 2)^2;
A = V * SIGMA * V';
b = A * xopt;
endLinear_solu.m
function [epoch, x_now] = linear_solu(A, b, x0, epsilon)
x_now = x0; epoch = 0;
gnorm = norm(A * x0 - b, 2);
ginit = gnorm;
while gnorm / ginit > epsilon
g_now = A * x_now - b;
gnorm = norm(g_now, 2);
xnext = x_now - (gnorm^2) / (g_now' * A * g_now) * g_now;
x_now = xnext;
epoch = epoch + 1;
end;
endLinear_impr.m
function [epoch, x_now] = linear_impr(A, b, x0, epsilon)
lr = @(epoch) - epoch / 5e4 + 0.9;
x_now = x0; epoch = 0;
gnorm = norm(A * x0 - b, 2);
ginit = gnorm;
while gnorm / ginit > epsilon
g_now = A * x_now - b;
gnorm = norm(g_now, 2);
xnext = x_now - lr(epoch) * (gnorm^2) / (g_now' * A * g_now) * g_now;
x_now = xnext;
epoch = epoch + 1;
end;
endLinear_eval.m
function linear_eval(type)
% linear_eval('normal') for normal gradient descent algorithm
% linear_eval('improv') for improved radient descent algorithm
epsilon = 1e-6;
fprintf('%8s ', 'EPOCH/ERR');
for n = 1: 5
fprintf(' n=%3d ', 100 * n);
end;
fprintf('n');
for c = 3: 6
fprintf('COND=1e+%d ', c);
for n = 1: 5
[A, b, xopt, x0] = linear_init(10^c, 100 * n);
if strcmp(type, 'improv')
[epoch, xk] = linear_impr(A, b, x0, epsilon);
else
[epoch, xk] = linear_solu(A, b, x0, epsilon);
end;
error = norm(xk - xopt, 2) / norm(x0 - xopt, 2);
fprintf('%5d/%.4f ' , epoch, error);
end;
fprintf('n');
end;
end实验结果
1
EPOCH/ERR n= 100 n= 200 n= 300 n= 400 n= 500
COND=1e+3 3382/0.0786 8910/0.0682 14400/0.0511 15750/0.0492 15696/0.0525
COND=1e+4 3382/0.0786 8910/0.0682 14400/0.0511 15750/0.0492 15696/0.0525
COND=1e+5 3382/0.0786 8910/0.0682 14400/0.0511 15750/0.0492 15696/0.0525
COND=1e+6 3382/0.0786 8910/0.0682 14400/0.0511 15750/0.0492 15696/0.0525
-epoch / 5e4 + 0.9
EPOCH/ERR n= 100 n= 200 n= 300 n= 400 n= 500
COND=1e+3 203/0.0786 459/0.0682 419/0.0510 588/0.0494 544/0.0536
COND=1e+4 221/0.0786 534/0.0682 855/0.0513 538/0.0494 761/0.0535
COND=1e+5 203/0.0786 333/0.0682 617/0.0509 538/0.0494 761/0.0535
COND=1e+6 203/0.0786 534/0.0682 617/0.0509 588/0.0494 544/0.0536
np.exp(- np.log(0.9) * (epoch / 4e3 - 1))
EPOCH/ERR n= 100 n= 200 n= 300 n= 400 n= 500
COND=1e+3 287/0.0786 462/0.0682 513/0.0512 579/0.0489 545/0.0535
COND=1e+4 313/0.0786 531/0.0682 705/0.0512 617/0.0496 710/0.0528
COND=1e+5 287/0.0786 412/0.0682 619/0.0510 617/0.0496 710/0.0528
COND=1e+6 287/0.0786 531/0.0682 619/0.0510 579/0.0489 545/0.0535算法变体
这里清川对原始算法做了一些小改进，在
这个改进并不能让人喜悦，显然它是在当前数据（矩阵A与b的值）下过拟合的。对于新的数据，该算法未必能够加速。但是这说明了另一个问题：最速下降法并不是最速的。
这是因为这里选取的欧式范数（
5956
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
