精华内容
下载资源
问答
  • 周期函数迭代法
    千次阅读
    2019-01-12 19:25:40

    梯度下降法、牛顿法、拟牛顿法 三类迭代法应用场景有何差别?

    By Datawhale知乎内容输出小组D1

    问题

    梯度下降法一族(如SGD、Adam)、牛顿法一族(如Gauss-Newton Method,LM法)、拟牛顿法一族(如L-BFGS)是机器学习中最常见的三大类迭代法,但三者分别通常擅长解决的应用场景是什么?为什么会这样的呢?谢谢

    解答

    梯度下降法(SGD为例)牛顿法拟牛顿法
    时间复杂度(单次迭代)只需计算1阶导,时间复杂度低,为 O ( n ) O\left( n \right) O(n)需计算Hessian矩阵及其逆,时间复杂度高,为 O ( n 3 ) O\left( {{n^3}} \right) O(n3)用正定矩阵近似Hessian矩阵的逆,时间复杂度为 O ( n 2 ) O\left( {{n^2}} \right) O(n2)
    收敛速度收敛慢,迭代次数大收敛快,迭代次数小-
    初始值要求无太强要求,容易逃离鞍点对初始值有一定要求,非凸问题容易陷入鞍点(牛顿法步长会越来越小)-
    应用场景特征维度较大的场景,如特征数>10k特征维度较小的场景需满足拟牛顿条件,更适合凸问题

    此外,在神经网络(非凸问题)的训练中,大多数都采用梯度下降法一族方法。而在训练逻辑回归(凸问题)等模型时,可采用梯度下降和拟牛顿方法。

    关于时间复杂度和收敛速度的差异,起因于求解方法:

    机器学习的任务中,是要最小化损失函数 L ( θ ) L\left( \theta \right) L(θ),其中 θ \theta θ 是待求的模型参数。梯度下降法、牛顿法/拟牛顿法都是迭代求解。梯度下降法是梯度求解,而牛顿法/拟牛顿法是用二阶的海森矩阵的逆矩阵或伪逆矩阵求解。

    迭代公式 θ t = θ t − 1 + Δ θ {\theta ^t} = {\theta ^{t - 1}} + \Delta \theta θt=θt1+Δθ

    求解方法
    梯度下降法:一阶泰勒展开
    L ( θ t ) = L ( θ t − 1 + Δ θ ) ≈ L ( θ t − 1 ) + L ′ ( θ t − 1 ) Δ θ \begin{aligned} L\left( {{\theta ^t}} \right) &= L\left( {{\theta ^{t - 1}} + \Delta \theta } \right) \\ &\approx L\left( {{\theta ^{t - 1}}} \right) + L'\left( {{\theta ^{t - 1}}} \right)\Delta \theta \end{aligned} L(θt)=L(θt1+Δθ)L(θt1)+L(θt1)Δθ

    为使 L ( θ t ) &lt; L ( θ t − 1 ) L\left( {{\theta ^t}} \right) &lt; L\left( {{\theta ^{t - 1}}} \right) L(θt)<L(θt1),可取 Δ θ = − α L ′ ( θ t − 1 ) \Delta \theta = - \alpha L&#x27;\left( {{\theta ^{t - 1}}} \right) Δθ=αL(θt1),则 θ t = θ t − 1 − α L ′ ( θ t − 1 ) {\theta ^t} = {\theta ^{t - 1}} - \alpha L&#x27;\left( {{\theta ^{t - 1}}} \right) θt=θt1αL(θt1)
    其中 α \alpha α 是步长,一般直接赋一个较小的值。

    牛顿法:二阶泰勒展开
    L ( θ t ) = L ( θ t − 1 + Δ θ ) ≈ L ( θ t − 1 ) + L ′ ( θ t − 1 ) Δ θ + L ′ ′ ( θ t − 1 ) Δ θ 2 2 \begin{aligned} L\left( {{\theta ^t}} \right) &amp;= L\left( {{\theta ^{t - 1}} + \Delta \theta } \right) \\ &amp;\approx L\left( {{\theta ^{t - 1}}} \right) + L&#x27;\left( {{\theta ^{t - 1}}} \right)\Delta \theta + L&#x27;&#x27;\left( {{\theta ^{t - 1}}} \right){{\Delta {\theta ^2}} \over 2} \end{aligned} L(θt)=L(θt1+Δθ)L(θt1)+L(θt1)Δθ+L(θt1)2Δθ2

    为简化分析过程,假定 θ \theta θ 只有一维,将一阶和二阶导数分别记为 g g g h h h,即:
    L ( θ t ) ≈ L ( θ t − 1 ) + g Δ θ + h Δ θ 2 2 L\left( {{\theta ^t}} \right) \approx L\left( {{\theta ^{t - 1}}} \right) + g\Delta \theta + h{{\Delta {\theta ^2}} \over 2} L(θt)L(θt1)+gΔθ+h2Δθ2
    为使 L ( θ t ) L\left( {{\theta ^t}} \right) L(θt) 极小,即 g Δ θ + h Δ θ 2 2 g\Delta \theta + h{{\Delta {\theta ^2}} \over 2} gΔθ+h2Δθ2 极小,令
    ∂ ( g Δ θ + h Δ θ 2 2 ) ∂ ( Δ θ ) = 0 ⇒ Δ θ = − g h {{\partial \left( {g\Delta \theta + h{{\Delta {\theta ^2}} \over 2}} \right)} \over {\partial \left( {\Delta \theta } \right)}} = 0 \Rightarrow \Delta \theta = - {g \over h} (Δθ)(gΔθ+h2Δθ2)=0Δθ=hg


    θ t = θ t − 1 − g h {\theta ^t} = {\theta ^{t - 1}} - {g \over h} θt=θt1hg

    θ \theta θ 推广到向量形式,迭代公式为
    θ t = θ t − 1 − H − 1 g {\theta ^t} = {\theta ^{t - 1}} - {H^{ - 1}}g θt=θt1H1g

    此处 H H H 为海森矩阵。

    Reference

    1.梯度下降、牛顿法、拟牛顿法
    2.梯度下降法与牛顿法比较

    更多相关内容
  • 数学实验第二课:函数迭代

    千次阅读 2020-07-26 15:39:04
    迭代 定理 在迭代函数连续的条件下,如果迭代序列...用迭代法近似计算 计算三次根号2 过程分析 首先取 f(x)=2/x^2 因为f(2(1/3))=2(1/3) 在该点的导数=2>1,不满足要求 >> syms x; >> y=diff(2/x^2);

    迭代

    在这里插入图片描述
    定理
    在这里插入图片描述
    在迭代函数连续的条件下,如果迭代序列收敛,则一定收敛于方程x=f(x)的根x*,这个根为函数的不动点。
    不动点处的一阶导数小于1,则称不动点x为稳定的,若一阶导数等于0,则称不动点为超稳定,在超稳定附近,迭代过程x(n+1)=f(x(n))收敛到x*是非常快的。

    用迭代法近似计算

    计算三次根号2
    过程分析
    首先取
    f(x)=2/x^2
    因为f(2(1/3))=2(1/3)
    在该点的导数=2>1,不满足要求

    >> syms x;
    >> y=diff(2/x^2);
    >> subs(y,x,2^(1/3))%用值代替x
    >> z=subs(y,x,2^(1/3));
    >> fprintf('z = %.2f\n',z)%保留两位小数
    z = -2.00
    

    对上式恒等变形
    x3=2
    x3+ax2+bx=ax2+bx+2
    x=(ax2+bx+2)/(x2+ax+b)
    取f(x)=(ax2+bx+2)/(x2+ax+b)
    令a=0,b=1,得f(x)=(x+2)/(x2+1)

    对新的方程求解

    >> syms x;
    >> y=(x+2)/(x^2+1);
    >> p=diff(y);
    >> subs(p,x,2^(1/3))
    >> z=subs(p,x,2^(1/3))
    >> fprintf('z = %.2f\n',z)
    z = -0.84
    

    但这样收敛速度有点慢

    继续变形
    令a=0,b=4,得f(x)=(4x+2)/(x2+4),
    f(x)=(ax+b)/(cx+d)
    要满足a=d,b=3c

    >> y=(4*x+2)/(x^2+4);
    p=diff(y);
    z=subs(p,x,2^(1/3)); fprintf('z = %.2f\n',z)
    z = 0.15
    
    >> f=inline('(4*x+2)/(x^2+4)');
    >> x0=12;%初值随便给一个
    >> for i =1:100
    x0=f(x0);
    fprintf('%g,%g\n',i,x0);
    end
    1,0.337838
    2,0.814595
    3,1.12754
    4,1.23501
    5,1.25606
    6,1.25935
    7,1.25984
    8,1.25991
    9,1.25992
    10,1.25992
    11,1.25992
    12,1.25992
    13,1.25992
    14,1.25992
    

    迭代上式,稳定在1.25992,即为其方程的解,即三次根号2的值为1.25992
    在这里插入图片描述
    验算logistics的收敛性
    例如f(x)=ax(1-x),(0<=x<=1),a=3.2,x0=0.2
    先编写一个Demo函数

    function Demo(f,x0,N)
    close all
    for i =1:N
        axis([50,N,0,1]);
        if i>50
            plot(i,f(x0),'.');
            hold on ;pause(0.05);%画图加延迟
        end 
        x0=f(x0);
    end;
    hold off
    clear all;
    

    命令行下敲入

    Demo(@(x)3.2*x*(1-x),0.2,200)
    

    在这里插入图片描述
    可以看出形成了循环,值一直在两个值上跳动,周期为2

    更改a的值,可以出现不同的收敛情况
    a=3.6,x0=0.6
    在这里插入图片描述
    出现混沌现象

    logistics模型

    f(x)=ax(1-x) (0<=x<=1)
    a>0,为logistics参数

    在logistic函数中,首先取a的值为3,在(0,1)中随机取一数x。作为初值,进行迭代,共迭代300次左右,丢弃起始的100次迭代的数据,在图上绘出所有的点(a,Xn)(>100).然后慢慢地增加a值,每增加一次,都重复前面的步骤,一直增加到a=4为止,画出最后所得到的图形。

    logistic=inline('u*x*(1-x)');
    for u=3.0:0.01:4
    x0=rand;
    for i =1:300
    x0=logistic(u,x0);
    if i>100
    plot(u,x0,'k','linewidth',1);
    hold on;
    end;
    end;
    end;
    hold off
    

    在这里插入图片描述
    形成了周期的循环
    画平行于y轴的直线

    在这里插入图片描述
    初值影响很大,所以叫做蝴蝶效应

    在这里插入图片描述

    close all;clear all;
    f=inline('3.99*x*(1-x)');
    N=300;x0_1=0.663489000;x0_2=0.663489001;
    for i=1:N
    axis([1,N,0,1]);
    yn=abs(f(x0_1)-f(x0_2));%取两个初值函数差的绝对值
    plot(i,yn,'.');
    hold on;pause(0.005);
    x0_1=f(x0_1);
    x0_2=f(x0_2);
    end;
    hold off;
    

    迭代初期,差距不大,但是后面产生混沌现象,导致结果差距很大,蝴蝶效应

    在这里插入图片描述
    周期为2,可以取Feigenbaum上的值找周期

    二维迭代

    f=@(x,y)y-sin(x);
    g=@(x)3.1-x;
    xn=1.2;
    yn=0;
    for n=1:1000
    xN=xn;
    yN=yn;
    xn=f(xN,yN);
    yn=g(xN);
    plot(xn,yn,'k*');
    axis([-4,6,-3,7]);
    pause(0.1);
    hold on;
    end
    hold off
    

    在这里插入图片描述
    小图形经过旋转、折叠后可以重合,这类图形是分形图

    Martin迭代

    function Martin(a,b,c,N)%N为迭代次数
    f=@(x,y)(y-sign(x)*sqrt(abs(b*x-c)));%f和g称为Martin迭代
    g=@(x)(a-x);
    m=[0;0];
    for n=1:N
        m(:,n+1)=[f(m(1,n),m(2,n)),g(m(1,n))];
    end
    plot(m(1,:),m(2,:),'kx');
    axis equal
    

    在这里插入图片描述
    局部与整体具有某种相似
    但不能从二维迭代散点图中判断迭代序列是否收敛

    展开全文
  • 方程组求解的直接法与迭代法实现

    千次阅读 2018-12-14 20:25:01
    方程组求解的直接法与迭代法实现 问题描述 我们的目的在于求解如下所示的方程组: 其中的A11、A12、A21、A22A_{11}、A_{12}、A_{21}、A_{22}A11​、A12​、A21​、A22​以及右端项fff以不同的的稀疏存储方式存在...

    方程组求解的直接法与迭代法实现

    问题描述

    我们的目的在于求解如下所示的方程组:

    在这里插入图片描述

    其中的 A 11 、 A 12 、 A 21 、 A 22 A_{11}、A_{12}、A_{21}、A_{22} A11A12A21A22以及右端项 f f f以不同的的稀疏存储方式存在文件中。

    我们需要完成的问题包括:

    • 理解并描述 A 11 、 A 12 、 A 21 、 A 22 A_{11}、A_{12}、A_{21}、A_{22} A11A12A21A22以及右端项 f f f的数据存储方式,其中 A 11 A_{11} A11是以CSR(compressed)方式存储,其余以COO(coordinate or IJ)方式存储。
    • 设计一种数据存储方式来储存具有块结构的稀疏矩阵 A A A
    • 写一段程序读取 A A A的子块,并将其以设计好的结构存储为 A A A
    • 用MUMPS,SuperLU等方法来直接求解方程组。
    • 选择和应用一种迭代方法来求解线性方程组,并提供随着迭代步的进行的收敛情况(和直接法做比较)。
    • 基于你对这个线性系统的理解,设计一种有效的预条件子方法来加速迭代方法的收敛。

    CSR和COO矩阵压缩稀疏存储方式说明

    CSR:压缩稀疏行存储

    N = 58324,M=99。
    A 11 A_{11} A11稀疏存储文件说明:

    • 第1行3个数:
      行块数(58324),列块数(58324),非零块数(385766)
    • 第2行:
      每一块小矩阵的规模(边长3)
    • 第3行:
      每一块的flatten方式(0表示一行一行来)
    • 第4行:
      记录非零矩阵块累计数(按行)的向量的元素个数(58325)
    • 第5~58329行(58325个):
      表示非零矩阵块数目的累计值(IA),即后一个值减前一个值表示当前行的非零块个数。
    • 第58330行:
      非零块的个数(385766)
    • 第58331~444096行(385766个):
      表示对应于(IA)的非零块的列标(JA,标号从0开始)
    • 第444097行:
      非零元素块的元素个数综合(3471894)
    • 第444098~3915991行(3471894个):
      非零元素位置对应的值
    COO:坐标稀疏存储

    A 12 A_{12} A12稀疏存储文件说明:

    • 第1行3个数:
      行数(174972),列数(99),非零元素(1148)
    • 第2到1149行(1148个):
      行坐标,列坐标,值

    A 21 A_{21} A21稀疏存储文件说明:

    • 第1行3个数:
      行数(99),列数(174972),非零元素(1100)
    • 第2到1101行(1100个):
      行坐标,列坐标,值

    A 21 A_{21} A21稀疏存储文件说明:

    • 第1行3个数:
      行数(99),列数(99),非零元素(99)
    • 第2到100行(99个):
      行坐标,列坐标,值

    f f f稀疏存储文件说明:

    • 第1行:
      列数(175071)
    • 第2到175072行(175071个):
      位置和值

    A A A的组装和存储

    观察 A A A的四个子块以及右端项的存储方式,我们发现,除了 A 11 A_{11} A11使用CSR方式存储之外,其他部分用的都是COO的存储方式。一个简单的想法是,先把 A 11 A_{11} A11装成COO存储方式,在将几个子块读入,组装成 A A A,以COO的方式写入文件,以便备用。这个过程我用Matlab(Matlab 2014a 盗版)来实现。

    A 11 A_{11} A11的从CSR格式转成COO格式并保存的代码如下:

    clc
    clear
    %% 读入数据,只读第一列
    filename = 'A11.dat';
    delimiter = ' ';
    formatSpec = '%f%*s%*s%[^\n\r]';
    fileID = fopen(filename,'r');
    dataArray = textscan(fileID, formatSpec, 'Delimiter', delimiter, 'MultipleDelimsAsOne', true, 'EmptyValue' ,NaN, 'ReturnOnError', false);
    fclose(fileID);
    A11 = [dataArray{1:end-1}];
    clearvars filename delimiter formatSpec fileID dataArray ans;
    %% 处理数据,转为COO格式,并保存
    n = 174972;
    IA = A11(5:58329,:);
    nonzero_mat_num = IA(2:end)-IA(1:end-1);%共58324个
    JA = A11(58331:444096,:);
    JA = JA+1;%matlab的下标从1开始
    values = A11(444098:3915991,:);
    row1 = [174972 174972 3471894];
    submat_sparse_store_row = [];
    for i=1:58324
        current_num = nonzero_mat_num(i);
        submat_sparse_store_row(end+1:end+current_num,1) = repmat(i,current_num,1);
    end
    submat_sparse_store_col = JA;
    sparse_X = [];
    sparse_Y = [];
    sparse_values = values;
    for i=1:385766
        ind_row = submat_sparse_store_row(i);
        ind_col = submat_sparse_store_col(i);
        [x,y] = meshgrid((ind_row-1)*3+1:(ind_row)*3,(ind_col-1)*3+1:(ind_col)*3);
        x_temp = x;
        y_temp = y;
        x_flatten = x_temp(:);
        y_flatten = y_temp(:);
        sparse_X(end+1:end+9,1) = x_flatten;
        sparse_Y(end+1:end+9,1) = y_flatten;
    end
    A11_sparse = sparse(sparse_X,sparse_Y,values);
    A11_sparse(2000,2000)
    sum(sum((A11_sparse ~= 0)))
    data = [row1;sparse_X-1 sparse_Y-1 values];
    save('A11_coo.dat','data','-ascii');%ps:save方式的存储会产生误差,但是速度快咯,anyway,问题不大
    %dlmwrite('A11_coo.dat',data,'delimiter','\t','precision','%0.6e');
    
    

    从四个子块组装为 A A A并以COO格式保存为A.dat的代码如下:

    % 将四个系数矩阵拼成一个矩阵
    clc
    clear
    %%
    N = 58324;
    M = 99;
    n = 3*N;
    m = M;
    %% 
    load('A11_coo.dat');
    load('A12.dat');
    load('A21.dat');
    load('A22.dat');
    
    %%
    A11_h = A11_coo(1,:);
    T11 = A11_coo(2:end,:);
    %%
    A12_h = A12(1,:);
    T12 = A12(2:end,:);
    T12(:,2) = T12(:,2)+n;
    %%
    A21_h = A21(1,:);
    T21 = A21(2:end,:);
    T21(:,1) = T21(:,1)+n;
    %%
    A22_h = A22(1,:);
    T22 = A22(2:end,:);
    T22(:,1) = T22(:,1) + n;
    T22(:,2) = T22(:,2) + n;
    
    %%
    A_header_num = A11_h(3) + A12_h(3) + A21_h(3) + A22_h(3);
    A_header = [175071 175071 A_header_num];
    T = [T11;T12;T21;T22];
    A = [A_header;T];
    %save('A.dat','A',''-ascii','-double');
    
    %% 写入数据,第一行和第一列和第二列以int类型写入,第三列(除第一行),以双精度类型写入。
    fileID = fopen('A.dat','w');
    fprintf(fileID,'%6d  %6d  %6d\r\n',A(1,1),A(1,2),A(1,3));
    %fprintf('%6d  %6d  %6d\n',A(1,1),A(1,2),A(1,3));
    for k=1:size(A,1)-1
    fprintf(fileID,'%d  %d  %1.6e\n',A(k+1,1),A(k+1,2),A(k+1,3));
    %fprintf('%d  %d  %1.6e\n',A(k+1,1),A(k+1,2),A(k+1,3));
    
    end
    fclose(fileID);
    
    
    %%
    % B = sparse(A(2:end,1)+1,A(2:end,2)+1,A(2:end,3));
    % B_full = full(B);
    

    Matlab的back slash方法求解方程组

    \符号求解方程组

    既然用了Matlab,首先,一个自然的想法就是使用Matlab的\来求解方程组。使用的代码如下:

    clc
    clear 
    
    %% 导入A和f的数据
    load('A.dat');
    f_id = fopen('f.dat','r');
    fgets(f_id);
    formatSpec = '%d %f';
    sizef = [2 Inf];
    ff = fscanf(f_id,formatSpec,sizef);
    fclose(f_id);
    ff = ff';
    f = ff(:,2);
    %%
    A_h = A(1,:);
    T = A(2:end,:);
    %%
    AA = sparse(T(1:end,1)+1,T(1:end,2)+1,T(1:end,3));
    u = AA\f;
    %%
    u_header = 175071;
    u_mine = u;
    save('u_mine.dat','u_mine','-ascii');
    
    %% check result
    E = u(1)-1.24456;
    
    
    

    做个简单的计时,发现back slash求解大概需要2.655s。

    在这里插入图片描述

    求解的结果的前几个值如下:

    在这里插入图片描述

    我的计算机参数

    自报一下家门,老电脑了,四核处理器,每个核心2.60GHz,单核每个时钟周期执行的浮点运算次数大约为25。
    理 论 峰 值 速 度 = cpu 主 频 ∗ 每 个 时 钟 周 期 执 行 的 浮 点 运 算 次 数 ∗ 核 数 理论峰值速度=\text{cpu}主频*每个时钟周期执行的浮点运算次数*核数 =cpu
    笔记本的峰值计算能力为256亿浮点运算/秒,具体如下:

    制造商 : Intel
    型号 : Intel® Core™ i5-3230M CPU @ 2.60GHz
    速度 : 3.2千兆赫(GHz)
    最低/最高速度/涡轮速度 : 1.2千兆赫(GHz) - 2.6千兆赫(GHz) - 3.2千兆赫(GHz)
    恒速 : 2.6千兆赫(GHz)
    峰值处理性能 (PPP) : 25.6数十亿浮点运算/秒(GFLOPS)
    调整后的峰值性能 (APP) : 7.68WG
    内核/处理器 : 2 个
    线程/内核 : 2 个
    类型 : 便携电脑

    关于matlab内置算法的一个说明

    关于matlab的backslash:对于稀疏矩阵的\(mldivide)方法求解,Matlab大多时候使用的是分解算法,如QR solver、diagonal solver、banded solver、(permuted)triangular solver、LU solver、Cholesky solver和LDLsolver等等。Matlab会自己进行条件判断,而根据矩阵情况选择合适的算法。

    另外,查看了matlab的document,其对于稀疏矩阵的分解解法和迭代解法如下:

    在这里插入图片描述

    对于一般的方法,有一定的适用条件:比如说pcg用于处理Hermitian正定矩阵,minres一般用于处理Hermitian矩阵。

    直接法求解(MUMPS, SuperLU,etc)

    Matlab没有现成的MUMPS,SuperLU的求解器。怎么办?只能使用c中的相应的库来实现。当然,也可以利用Matlab和c或者Fortran的混编,毕竟MUMPS的SuperLU的包没有matlab版本的,那么又想用Matlab来调用,这么干也是可以的。

    我第一个想到的是使用所里的phg来实现,因为它有相关的接口。我在ubuntu下安装了phg及其相关的依赖,用以调试程序,等到合适的时候,再将程序放到所里的科学计算集群(LSSC-IV)上去运行。

    phg的方程组求解器接口调用

    编写c程序,读入Matlab生成的A.dat和已有的f.dat文件,调用phg的解法器接口,进行求解。代码如下:

    #include "phg.h"
    #include <string.h>
    #include <math.h>
    #include <stdio.h>
    #include <stdlib.h>
    
    int main(int argc, char *argv[])
    {
    
    
    	/* Read data, including A and f. 先一次性把数据读进来,这样比较节省时间。*/
    	FILE *fp;
    	if ((fp = fopen("A.dat", "r")) == NULL)
    	{
    		printf("cannot find what you need,please check it!");
    		exit(0);
    	}
    	int m, n, numOfNonzero;
    	int i;
    	fscanf(fp, "%d%d%d", &m, &n, &numOfNonzero);
    	float *valuesOfA = (float*)malloc(sizeof(float)*numOfNonzero);
    	int *rows = (int*)malloc(sizeof(int)*numOfNonzero);
    	int *cols = (int*)malloc(sizeof(int)*numOfNonzero);
    	for (i = 0; i <= numOfNonzero - 1; i++)
    	{
    		fscanf(fp, "%d%d%f", &rows[i], &cols[i], &valuesOfA[i]);
    	}
    	fclose(fp);
    
    
    	if ((fp = fopen("f.dat", "r")) == NULL)
    	{
    		printf("cannot find what you need,please check it!");
    		exit(0);
    	}
    	int numOff;
    	fscanf(fp, "%d", &numOff);
    	double *valuesOff = (double*)malloc(sizeof(double)*numOff);
    	int *indexOff = (int*)malloc(sizeof(int)*numOff);
    	for (i = 0; i <= numOff - 1; i++)
    	{
    		fscanf(fp, "%d%lf", &indexOff[i], &valuesOff[i]);
    		//printf("%d,%0.15e\n", indexOff[i], valuesOff[i]);
    	}
    	fclose(fp);
    	/*finish reading data. */
    	//	printf("%d,%d,%0.6e\n", rows[0], cols[0], valuesOfA[0]);
    	//	printf("%d,%0.15e\n", indexOff[0], valuesOff[0]);
    	//	getchar();
    
    //读如的数据存放在rows、cols、valuesOfA和indexOff以及valuesOff中。
    
    
    
    	SOLVER *solver;//创建解法器对象
    	MAT *A;
    	MAP *map;
    //	MPI_Comm comm;
    //	phgInitSetComm(comm);
    	phgInit(&argc,&argv);
    //	printf("phgComm=%d",phgComm);
    	map = phgMapCreateSimpleMap(phgComm,-1,175071);
    	A = phgMapCreateMat(map,map);
    	int numOfNonzeroInA = 3474241;
    
    	for (i = 0; i <= numOfNonzeroInA - 1; i++)
    	{
    		int row = rows[i];
    		int col = cols[i];
    		float value = (float)valuesOfA[i];
    //printf("row=%d,col=%d,value=%f\n",row,col,value);
    
    		phgMatAddGlobalEntry(A, row, col, value);
    //printf("debug: program comes here……\n");
    
    	}
    	
    	
    	solver = phgSolverMat2Solver(SOLVER_DEFAULT,A);
    	solver->verb = 1;
    
    	int N = 175071;
    	
    	//printf("debug: program comes here……%f%f",indexOff[100],valuesOff);
    //	phgVecAddGlobalEntries(u,0,N,indexOff,valuesOff);
    	
    	
    	phgVecDisassemble(solver->rhs);	
    	for (i = 0; i <= N - 1; i++)
    	{
    		int index = indexOff[i];
    		float value = valuesOff[i];
    		value = (float)value;
    	//	printf("i=%dindex=%dvalue=%f",i,index,value);
    	//	getchar();
    	//	phgMatAddGlobalEntry(A, row, col, value);
    		phgSolverAddGlobalRHSEntry(solver,index,value);
    
    	}
    
    	VEC *u_h;//定义VEC对象存放解
    	u_h = phgMapCreateVec(map,1);
    //	phgVecDisassemble(u_h);
    	
    //	const char *var_name = "f_out";
    //	const char *file_name = "f_out";
    //	phgVecDumpMATLAB(solver->rhs,var_name,file_name);
    
    //	var_name = "A_out";
    //	file_name = "A_out";
    //	phgMatDumpMATLAB(solver->mat,var_name,file_name);
    	
    	phgSolverVecSolve(solver, 0, u_h);//求解
    /*	const char *var_name = "u_h";
    	const char *file_name = "u_mine";
    	phgVecDumpMATLAB(u_h,var_name,file_name);
    */
    	return 0;
    	
    }	
    
     	/*自由度对象赋值*/
      /*  	//solver = phgSolverCreate(SOLVER_DEFAULT, u_h, NULL);//定义解法器
    	for (int i = 0; i <= numOfNonzero - 1; i++)
    	{
    		int row = rows[i];
    		int col = cols[i];
    		float value = (float)valuesOfA[i];
    		phgSolverAddGlobalMatrixEntry(solver, row, col, value);
    
    	}
    	printf("debug: program comes here……\n");
    
    	for (int i = 0; i <= numOfNonzero - 1; i++)
    	{
    		int index = indexOff[i];
    		float value = (float)valuesOff[i];
    		phgSolverAddGlobalRHSEntry(solver, index, value);
    	
    
    	}
    
     	phgSolverSolve(solver, TRUE, u_h, NULL);//求解
    	//SIMPLEX *e;
    	//	const FLOAT lambda[];
    	//	FLOAT *values;
    	//	phgDofEval(u_h, e, lambda, values);//不足的输入要补全。
    	//	phgPrintf("first value of u is : %lf\n", *values);//打印结果
    	phgSolverDestroy(&solver);//销毁解法器对象
    */
    
    
    
    

    编译生成可执行文件后,我们来命令行传入参数来调用phg中的解法器来求解这个 方程组。phg中预装的解法器接口有: “pcg”, “bicgstab”, “gmres”, “lgmres”, “preonly”, “aps”, “2grid”, “ams”, “gather”, “asm”, “smoother”, “hypre”, “petsc”, “mumps”, “superlu”, “minres”, "diagonal"等。

    按照要求,我现在想调用mumps和superlu求解。

    基于phg平台的mumps方法求解

    phg提供的mumps的接口,提供了以下的参数:

    在这里插入图片描述

    在矩阵的性质中,我们给它指定非对称,其他的参数保持默认。使用如下所示命令行参数执行可执行文件:

    ./equationSolverTemplate -solver mumps -mumps_symmetry unsym
    

    结果如下:

    • 分析步:
      在这里插入图片描述

    • 分解步:
      在这里插入图片描述

    • 求解和检查
      在这里插入图片描述

    可以看到,打印出来的墙上时间为1.66s,比Matlab略快一些。

    基于phg平台的superlu方法求解

    查看superlu的参数帮助信息如下:

    在这里插入图片描述

    我选择保持默认。使用如下命令执行可执行程序。

    ./equationSolverTemplate -solver superlu
    

    结果如下:

    在这里插入图片描述

    消耗的墙上时间约为3.54s,比MUMPS方法略慢一些。

    Anyway,我们由直接法,不管是Matlab的back slash还是MUMPS抑或是SUPER LU,我们用直接法得到了一个可以在迭代法中作为标准的近似真解。我们可以把这个解保存下来,下面再使用这个解来作为迭代法收敛性分析的一个标准。

    迭代方法设计与求解

    A A A矩阵性态的观察

    使用Matlab的spy函数,我们来观察一下矩阵 A A A的分布,结果如下:

    在这里插入图片描述

    看起来是个对角矩阵,而且很稀疏。我们再通过sum(sum((A - A')))判断,发现 A A A并不是对称的,只不过是在外形上看起来有些对称(事实上,在外形上也不是对称的)。进一步,我们通过[r,p]=chol(A)判定,它并不是正定的(非对称意义下的正定)。

    所以,这是一个非零元素基本上都集中在对角线以及两条边上的没有特别好的性质的稀疏矩阵。

    方法的选择

    我们可以设计和使用多层迭代方法来求解。时间的原因,我们依然调用别人已经写好的包,而不自己动手完成这个过程。能使用的包就比较多了,比如由张晨松老师等人编写维护的FASP的包,下载地址,另外还有Matlab内置的共轭梯度,残差法等等。

    因为这个矩阵的在性质上,好像除了看起来有点对称以及稀疏之外,没有太好的性质。所以,我想采用GMRES方法求解。GMRES中文叫广义极小剩余法,它是上世纪末到本世界初被人们认为是求解大型稀疏非对称线性方程组的最有效的方法之一。Matlab有gmres内置函数,基于重启动的gmres方法。我想用的依然是phg的内置的gmres算法。

    为了加快的速度,我用的是预条件子的GMRES方法。phg中gmres solver的参数如下所示。

    在这里插入图片描述

    预条件子的选择,使用GMRES,我的预条件子选择了两个:一个是MUMPS作为预条件子,另一个是ASM(加性 Schwarz作为预条件子),并且在ASM预条件子中使用单精度的稀疏直接求解器MUMPS求解子区域问题。

    基于MUMPS预优的GMRES方法

    指定用单精度的MUMPS解法器做为GMRES的预条件子。我在 Makefile中定义的命令行参数如下:

    在这里插入图片描述

    上图中,gmres_pc_type 中的 “solver” 参数表示在 GMRES 每步迭代中调用另一个解法器做为预条
    件子,而 -gmres_pc_opts 选项则用来设定预条件子解法器及参数。

    运行的结果如下:

    在这里插入图片描述

    使用单个进程,设定atol=0, rtol=1e-10, btol=1e-10。从结果中可以看出,迭代2次之后,有一次重启动。迭代4次残差相对值达到了-18次方量级的。墙上时间消耗时间为2.04s。收敛速度是非常快的。这是在最外层的迭代,内层表现我们就不往里看了。

    基于ASM预优(ASM中使用MUMPS求解子区域问题)的GMRES方法

    指定使用 GMRES 迭代、ASM (加性 Schwarz)预条件子,并且在 ASM 预条件子中使用单精度的稀疏直接求解器 MUMPS 求解子区域问题。

    在Makefile中定义的参数列表如下:

    在这里插入图片描述

    结果如下:

    在这里插入图片描述

    经过一次的迭代(最外层的求解器),误差(相对残差)达到了e-14量级的。墙上时间消耗为1.81s。

    使用Matlab的gmres内置函数求解

    我们考虑一种比较简单的情况,即使用不完全LU预优分解,不采用重启动技术,我们来看看它的收敛表现。

    这一次,我们不用残差来衡量误差,而使用直接法得到的解作为一个参考标准。写一段简单的matlab脚本如下,用来调用matlab的gmres内置函数。

    clc
    clear
    %% 导入A,f和真解u
    load('A.dat');
    %A_h = A(1,:);
    T = A(2:end,:);
    AA = sparse(T(1:end,1)+1,T(1:end,2)+1,T(1:end,3));
    A = AA;
    
    f_id = fopen('f.dat','r');
    fgets(f_id);
    formatSpec = '%d %f';
    sizef = [2 Inf];
    ff = fscanf(f_id,formatSpec,sizef);
    fclose(f_id);
    ff = ff';
    f = ff(:,2);
    
    load u_mine.dat
    u = u_mine;
    
    clear AA f_id ff formatSpec sizef T u_id uu;
    %%
    [L,U] = ilu(A,struct('type','ilutp','droptol',1e-6));
    l = length(f);
    tol = 1e-12;
    maxit = 20;
    b = f;
    
    x0 = zeros(l,1);
    errors = [norm(x0 - u)/norm(u);];
    for i = 1:maxit
    [x0,fl0,rr0,it0,rv0] = gmres(A,b,1,tol,1,L,U,x0);
    errors(end+1) = norm(x0 - u)/norm(u);
    end
    st = 4;
    semilogy(st:maxit,errors(st+1:end),'-o');
    xlabel('Iteration number');
    ylabel('Relative Error');
    

    结果如下:

    在这里插入图片描述

    这里定义的误差为 Error = ∣ ∣ u − u h ∣ ∣ 2 / ∣ ∣ u ∣ ∣ 2 \text{Error} = ||u-u_h||_2/||u||_2 Error=uuh2/u2。这里的u,用的是直接法得出来的解。可以看得出来,使用不完全LU分解做预优的GMRES方法,只要两三步就达到比较高的精度了。也容易看出,迭代到第4步,误差就降低到了e-8量级的。另外,这个过程也不是单调的,在第5步到第6步,误差有一个非常微小的升高。

    展开全文
  • 将得到的数组x,y,构建一个机器学习模型,采用梯度下降,通过多次迭代,学习到函数的系数。使用python和numpy进行编程,具体实现的代码如下:import numpy as np%matplotlib inlinefrom matplotlib import pyplot...

    给出一个数组x,然后基于一个二次函数,加上一些噪音数据得到另一组数据y。

    将得到的数组x,y,构建一个机器学习模型,采用梯度下降法,通过多次迭代,学习到函数的系数。使用python和numpy进行编程,具体实现的代码如下:

    import numpy as np

    %matplotlib inline

    from matplotlib import pyplot as plt

    np.random.seed(100)

    x=np.linspace(-1,1,100).reshape(100,1)

    y=3*np.power(x,2)+2+0.2*np.random.rand(x.size).reshape(100,1)

    plt.scatter(x,y)

    plt.show()

    w1=np.random.rand(1,1)

    b1=np.random.rand(1,1)

    lr=0.001

    for i in range(800):

    y_pred=np.power(x,2)*w1+b1

    loss=0.5*(y_pred-y)**2

    loss=loss.sum()

    grad_w=np.sum((y_pred-y)*np.power(x,2))

    grad_b=np.sum((y_pred-y))

    w1-=lr*grad_w

    b1-=lr*grad_b

    plt.plot(x,y_pred,'r-',label='predict')

    plt.scatter(x,y,color='blue',marker='o',label='true')

    plt.xlim(-1,1)

    plt.ylim(2,6)

    plt.legend()

    plt.show()

    print(w1,b1)

    原始数据如图所示:

    6cc7316d125f

    原始数据.png

    得到的结果如图所示:

    6cc7316d125f

    拟合结果.png

    补充:

    使用pytorch深度学习框架中的Antograd实现函数拟合

    import torch as t

    import torch

    %matplotlib inline

    from matplotlib import pyplot as plt

    t.manual_seed(100)

    dtype=t.float

    x=t.unsqueeze(torch.linspace(-1,1,100),dim=1)

    y=3*x.pow(2)+2+0.2*torch.rand(x.size())

    plt.scatter(x.numpy(),y.numpy())

    plt.show()

    w=t.randn(1,1,dtype=dtype,requires_grad=True)

    b=t.zeros(1,1,dtype=dtype,requires_grad=True)

    lr=0.001

    for ii in range(800):

    y_pred=x.pow(2).mm(w)+b

    loss=0.5*(y_pred-y)**2

    loss=loss.sum()

    loss.backward()

    with t.no_grad():

    w-=lr*w.grad

    b-=lr*b.grad

    w.grad.zero_()

    b.grad.zero_()

    plt.plot(x.numpy(),y_pred.detach().numpy(),'r-',label='predict')

    plt.scatter(x.numpy(),y.numpy(),color='blue',marker='o',label='true')

    plt.xlim(-1,1)

    plt.ylim(2,6)

    plt.legend()

    plt.show()

    print(w,b)

    原始数据画图:

    6cc7316d125f

    原始数据

    函数拟合结果:

    6cc7316d125f

    函数拟合结果

    展开全文
  • 输入函数和解出现的周期,然后第一次猜测和迭代次数
  • Jacobi(雅可比)迭代原理与matlab代码

    万次阅读 多人点赞 2020-10-26 21:57:58
    Jacobi(雅可比)迭代原理与matlab,C代码 Jacobi迭代分量形式:  xi(k+1)=1/aii(bi−∑j≠ii=1naijxj(k)) \ x_i^{(k+1)}=1/a_{ii}(b_i-\sum_{\stackrel{i=1}{ j\neq i} }^na_{ij}x_j^{(k)})  xi(k+1)​=1/...
  • 北航计算机学院-计算机组成原理课程设计-2020秋,PreProject-Logisim-Logisim组合逻辑电路-斐波那契数列问题(简单迭代法+矩阵乘法的快速幂)。 北航计算机学院的计算机组成原理课程设计,是高度实践性的专业课程,...
  • m文件“dpre”通过循环QZ或牛顿反向迭代法解决离散时间周期最优控制问题。 这些不是可用的最快方法,但效果很好。 mex 文件“dprex”通过周期性 QR(使用来自 matlab 内部 slicot 库的函数)或复杂的周期性 QC ...
  • 这样一来每个节点间的误差在加总时就不会相互抵消,但绝对值的存在使得函数图像会变成一个”V”字型,在最低点处是一个箭头,于是这个函数在最低点处不连续,梯度下降不能运用于不连续的函数。 第三者是两者...
  • 迭代最小二乘、牛顿及其matlab实现

    千次阅读 多人点赞 2021-04-10 20:19:37
    角度测量(AOA/DOA)室内定位-迭代最小二乘和高斯牛顿\MATLAB 原创不易,路过的各位大佬请点个赞 AOA定位,角度测量 迭代最小二乘、高斯牛顿 二维仅角度测量的定位问题; ~ ~ !!!!!!!!!!!!!!!!!...
  • 函数式编程中的重要概念

    千次阅读 2019-06-26 11:15:17
    函数式编程中的重要概念函数式编程范式的意义函数类型与高阶函数部分函数柯里化闭包递归记忆化 原文地址 函数式编程范式的意义 在众多的编程范式中,大多数开发人员比较熟悉的是面向对象编程范式。一方面是由于面向...
  • Python-函数基础总结与内置函数

    千次阅读 2021-01-13 15:57:47
    调用函数位置传参与关键字传参传参是值传递还是引用传递定义函数参数默认参数关键字参数参数组返回值指定参数、返回值类型内置函数标准类型函数dirhelpidlenstrtype数字类型函数转换工厂函数功能函数用于可迭代对象...
  • 与冲激函数、阶跃函数的卷积.ppt

    千次阅读 2021-04-22 01:49:33
    与冲激函数、阶跃函数的卷积信号与系统总 复 习 第一章 绪论 1、信号的概念 2、分类:典型的连续时间信号: 指数、正弦、复指数、抽样、钟形、δ(t), u(t), eat, sin(ω0t), Sa(kt) 3、信号的运算: 移位、反褶、...
  • 由于开普勒方程属于超越方程,因而无法通过解析精确求解,这一问题在历史上困扰全世界数学家们达 200 年之久,直至牛顿迭代方法的出现。本文将介绍一种实用的开普勒方程求解方法,并采用 C 语言实现其算法。该方法...
  • 在强化学习中,当环境模型已知时(也即环境状态转移概率和奖励已知),可以采用动态规划的思想来解决强化学习问题,常用的有策略迭代算法和值迭代算法两种,以下展开具体介绍。 1.动态规划与强化学习的联系 动态...
  • 目录概述神经网络感知器(Perceptron)多层感知器(MLP)常用函数 概述 在我们学习深度学习模型,或者了解某种算法的过程中,经常会看到如卷积层,全连接层,归一化,正则化,激活函数等概念,这些东西如果直接看的...
  • 作为一名MATLAB的狂热学习者,从最开始接触到至现在的日趋炉火纯青,确实是一段辛酸的成长史...diric 产生Dirichlet函数周期Sinc函数 gauspuls 产生高斯调制正弦脉冲 pulstran 产生脉冲串 rectpuls 产生非周期矩形...
  • 本发明属于数据分析技术领域,涉及筛选迭代余量的相对方差作为经验模态分解方法筛选迭代过程的终止准则。背景技术:一维的Fourier分解、小波分析,二维的PCA/EOF等方法,都是从低频开始分解,获取不同的模态/特征...
  • 而在C++中,初始化时在执行相关代码时才会进行初始化,主要是由于C++引入对象后,要进行初始化必须执行相应构造函数和析构函数,在构造函数或析构函数中经常会需要进行某些程序中需要进行的特定操作,并非简单地...
  • 本文就深度学习中经常使用的损失函数、优化器、学习率迭代策略等进行介绍。 一、损失函数 损失函数是指用于计算标签值和预测值之间差异的函数,对于不同的任务,使用的损失函数也不相同,在目标检测等领域,损失函数...
  • 角度测量(AOA/DOA)室内定位-迭代最小二乘和高斯牛顿\MATLAB 原创不易,路过的各位大佬请点个赞 二维仅角度测量的定位问题; 方法:1-非线性最小二乘估计-迭代最小二乘求解 2-极大似然估计-高斯牛顿求解 性能指标...
  • 深度学习基础宝典---激活函数、Batch Size、归一化

    千次阅读 多人点赞 2022-07-27 20:45:40
    深度学习基础宝典---激活函数、Batch Size、归一化
  • LASSO回归损失函数详解

    千次阅读 2020-03-11 14:37:48
    不过梯度下降和坐标轴下降的共性就都是迭代法,通过启发式的方式一步步迭代求解函数的最小值。 坐标轴下降法的数学依据主要是这个结论(此处不做证明):一个可微的凸函数J(θ), 其中θ是nx1的向量,即有n个维度。...
  • C语言自定义函数使用

    千次阅读 2021-08-26 23:05:56
    1.自定义函数(有函数名,返回值类型,函数参数) 2.函数的组成 ret_type fun_name(paral,*) ( statement; ) ret_type返回类型 fun_name函数名 paral 函数参数 int Add(int x,int y)//第一个int是返回类型 { //...
  • 迭代最小二乘和高斯牛顿\MATLAB Target localization with angle-only measurements 二维仅角度测量的定位问题; 方法:1-非线性最小二乘估计-迭代最小二乘求解 2-极大似然估计-高斯牛顿求解 性能指标:RMSE...
  • 直接利用 C 库函数中的标准整数除程序要花费 20~100 个周期,消耗较多资源。在非嵌入式领域,因为 CPU 运算速度快、存储器容量大,所以执行除运算和求模运算消耗的这些资源对计算机来说不算什么。但是在嵌...
  • Matlab优化函数中options选项的修改

    万次阅读 2018-04-25 23:43:13
    初学matlab优化,迭代中止后,经常一头雾水。参看帮助后仍似懂非懂。下面关于fminbnd函数的说明(也可作为fmincon函数的参考)对于新手也许会有帮助,不当之处请指正。 目标函数fun: 需要最小化的目标函数。fun函数...
  • 损失函数和梯度下降

    千次阅读 2021-03-12 11:03:10
    前向传播经过若干个神经元,再经过激活函数,最终得到结果,然后输出损失函数,根据损失函数再进行反向传播,及传递梯度,调整权重。 并不是根据激活偶函数输出直接返回梯度值,而是在计算损失函数的基础上进行反向...
  • GPS从入门到放弃(二十六) --- RTKLIB函数解析

    万次阅读 多人点赞 2020-06-07 21:57:21
    GPS从入门到放弃(二十六) — RTKLIB函数解析 为了贴合这个系列的标题“从入门到放弃”,在入门之后现在就要放弃此方向了。虽然感觉遗憾,暂时也没有办法。在此附上此系列最后一篇,希望能给大家一些帮助。 此文中...
  • Tensorflow常用函数汇总

    万次阅读 多人点赞 2018-03-01 08:54:02
    tf.mod(x, y, name=None) 取模 tf.abs(x, name=None) 求绝对值 tf.neg(x, name=None) 取负 (y = -x). tf.sign(x, name=None) 返回符号 y = sign(x) = -1 if x ; 0 if x == 0; 1 if x > 0. tf.inv(x, name=None)...

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 23,687
精华内容 9,474
热门标签
关键字:

周期函数迭代法

友情链接: PHPexcersize.rar