广州网站外包,网站专栏建设工作方案,黑河哈尔滨网站建设,清河做网站哪儿便宜凸二次优化问题Theory. 设是实对称正定矩阵#xff0c;#xff0c;则求解凸二次优化问题等价于求解线性方程组。Proof. 二次型二阶可导#xff0c;极小值点处梯度为零#xff0c;现对优化的目标函数求梯度。二次型本质上具有#xff1a;计算梯度的分量表达式#xff1a;合…凸二次优化问题Theory. 设是实对称正定矩阵则求解凸二次优化问题等价于求解线性方程组。Proof. 二次型二阶可导极小值点处梯度为零现对优化的目标函数求梯度。二次型本质上具有计算梯度的分量表达式合在一起写成矩阵形式显然凸二次优化问题的极值条件等价于该线性方程组。凸二次优化问题在建模中十分常见这说明讨论线性方程组的求解方法具有普遍的实用价值。然而对于规模较大的问题使用线性代数中的克莱姆法则暴力展开将导致时间开销巨大而高斯消元法算法流程又较为复杂。本文将介绍一种常见的数值分析方法求得线性方程组的数值近似解。最速梯度下降法称优化目标函数的梯度为残量Residue即是当前解对于线性方程组的不满足量由于函数是凸函数极值点一定存在。当前解处函数的梯度值表示了函数值上升最快的方向梯度方向上方向导数最大。那么沿着相反的方向每迭代一步就会更加靠近最优的极小值解。于是我们定义迭代关系其中表示迭代第轮的解表示每轮迭代的步长即每一步下降多少的权重。之所以需要设计一个与相关的步长是因为随着迭代的进行梯度是变化的。越靠近驻点的梯度将会越小直到接近于0时收敛。同时还要考虑在接近于驻点时应该放缓步伐否则会出现梯度在正负之间频繁震荡解在最优解左右摇摆的情况。在最速下降法中将作为自变量代入原函数并看作的函数对其进行最小化同样的驻点处梯度为零。根据链式求导法则最后一步由于是正定矩阵所以分母也是一个正定二次型值不为零可以直接除过来。由此我们最终得到了解的迭代关系算例设计与实验参考南科大的数值分析作业题EH计算设计SUSTech的算例一由于版权关系无传送门生成实对称正定矩阵 采用文献[1]中的类似方法生成矩阵其中和分别为Householder矩阵和对角矩阵Householder矩阵是对称正交矩阵这时对角矩阵的特征值就是的特征值参考矩阵的正交分解。取推导条件数计算公式 矩阵的最大特征值和最小特征值分别为那么其特征值介于之间对称正定矩阵的条件数公式为于是可以反解处条件数的计算公式为了能够重复试验数据使用线性同余法计算伪随机数进行向量的初始化也可以使用相同的随机种子调用编程语言内部的随机值库这里尊重南科大的原题要求。令xopt表示最优解x0表示解的初始向量。伪随机数算法参数给定如下面MATLAB程序所示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;这里设定最优解的分量在之间不为零的初始向量分量在之间。根据以上公式可以根据向量得到矩阵的值其后使用计算得到向量的值。这时将和作为问题的给定参数x0作为初始向量使用梯度下降的方式计算xk并计算它与最优解的真值之间的相对误差这里设定算法的停机标准为当前梯度变为接近于零的极小量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(fCOND1e{c}, end )for n in range(1, 6):A, b, xopt, x0 initialize(cond10**c, numb100*n)epoch, xk gradientDescent(A, b, x0, epsilon1e-6)error np.linalg.norm(xk - xopt) / np.linalg.norm(x0 - xopt)print({:5}/{:.4f}.format(epoch, error), end )print()MATLAB程序代码Linear_init.mfunction [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.mfunction [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.mfunction [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.mfunction 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(COND1e%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
COND1e3 3382/0.0786 8910/0.0682 14400/0.0511 15750/0.0492 15696/0.0525
COND1e4 3382/0.0786 8910/0.0682 14400/0.0511 15750/0.0492 15696/0.0525
COND1e5 3382/0.0786 8910/0.0682 14400/0.0511 15750/0.0492 15696/0.0525
COND1e6 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
COND1e3 203/0.0786 459/0.0682 419/0.0510 588/0.0494 544/0.0536
COND1e4 221/0.0786 534/0.0682 855/0.0513 538/0.0494 761/0.0535
COND1e5 203/0.0786 333/0.0682 617/0.0509 538/0.0494 761/0.0535
COND1e6 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
COND1e3 287/0.0786 462/0.0682 513/0.0512 579/0.0489 545/0.0535
COND1e4 313/0.0786 531/0.0682 705/0.0512 617/0.0496 710/0.0528
COND1e5 287/0.0786 412/0.0682 619/0.0510 617/0.0496 710/0.0528
COND1e6 287/0.0786 531/0.0682 619/0.0510 579/0.0489 545/0.0535算法变体这里清川对原始算法做了一些小改进在前面加上了一个人为设定的二次权重值。实验结果中每部分的第一行代表了的取值方式表中斜线前面代表停机时的总迭代次数后面代表相对误差。可以看出在保持解的精度不变的前提下加上人工修正的权重使得算法快了30倍。一开始清川想改进成随机梯度下降但显然他理解错了SGD算法的随机指的是神经网络训练时每轮迭代选取随机的样本和这里没有什么关系。但是清川从SGD中采纳了指数衰减和线性衰减的方法并加以调参就得到了改进。这个改进并不能让人喜悦显然它是在当前数据矩阵A与b的值下过拟合的。对于新的数据该算法未必能够加速。但是这说明了另一个问题最速下降法并不是最速的。这是因为这里选取的欧式范数分子上的二范数并不一定是最合适的衡量标准。每轮迭代总是假设梯度在当前迭代时是近似不变的以初始位置的梯度模拟整个步长上的梯度这样并不准确。而清川加了人工调参的二次权重后更好地模拟了梯度的变化所以更快。如果想进一步实验这一点可以将直接用替换去调参预计仍然能得到一个很好的收敛速度。做实验前最好对进行归一化否则调参过程中可能出现数值溢出。关于最速下降法不总是最速的这个结论也可以参考梯度下降法和最速下降法区别。注意本文并未对这两个概念作区分。
