凸二次优化问题
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 npdef 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 = 0for j in range(numb):t = (t * 31416 + 13846) % 46261v[j] = t * (2 / 46261) - 1t = 0for j in range(numb):t = (t * 42108 + 13846) % 46273xopt[j] = t * (5 / 46273) + 5for j in range(numb):t = (t * 42108 + 13846) % 46273x0[j] = t * (5 / 46273) - 10V = np.identity(numb) - 2 * v @ v.T / np.linalg.norm(v) ** 2A = V @ SIGMA @ V.Tb = A @ xoptreturn A, b, xopt, x0def gradientDescent(A, b, x0, epsilon):x_now = x0; epoch = 0gnorm = ginit = np.linalg.norm(A @ x0 - b)while gnorm / ginit > epsilon:g_now = A @ x_now - bgnorm = np.linalg.norm(g_now)xnext = x_now - (gnorm ** 2) / (g_now.T @ A @ g_now) * g_now # * LR_linearDecay(epoch)x_now = xnextepoch += 1return epoch, x_nowdef LR_linearDecay(epoch):return - epoch / 5e4 + 0.9def 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: numbt1 = mod(t1 * 31416 + 13846, 46261); v(j) = t1 * (2 / 46261) - 1; end; for j = 1: numbt2 = mod(t2 * 42108 + 13846, 46273); xopt(j) = t2 * (5 / 46273) + 5; end;for j = 1: numbt2 = 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 > epsilong_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 > epsilong_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 algorithmepsilon = 1e-6;fprintf('%8s ', 'EPOCH/ERR');for n = 1: 5fprintf(' n=%3d ', 100 * n);end;fprintf('n');for c = 3: 6fprintf('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.0536np.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的值）下过拟合的。对于新的数据，该算法未必能够加速。但是这说明了另一个问题：最速下降法并不是最速的。
这是因为这里选取的欧式范数（
...)
