精华内容
下载资源
问答
  • CG方法 CG标准函数库 万次阅读 多人点赞
    2016-10-13 20:37:21
     本文出处:http://blog.csdn.net/lcbcsdn/article/details/46848125
    

    (1)数学函数

    函数 功能描述
    abs(x) 返回输入参数的绝对值
    acos(x) 反余切函数,输入参数范围为[-1,1], 返回[0,π]区间的角度值
    all(x) 如果输入参数均不为0,则返回ture; 否则返回flase。&&运算
    any(x) 输入参数只要有其中一个不为0,则返回true。
    asin(x) 反正弦函数,输入参数取值区间为 [1,1] ,返回角度值范围为,  [π2,π2]
    atan(x) 反正切函数,返回角度值范围为 [π2,π2]
    atan2(y,x) 计算y/x的反正切值。实际上和atan(x)函数功能完全一样,至少输入参数不同。atan(x) = atan2(x, float(1))。
    ceil(x) 对输入参数向上取整。例如: ceil(float(1.3)) ,其返回值为2.0
    clamp(x,a,b) 如果x值小于a,则返回a
    如果x值大于b,返回b
    否则,返回x
    cos(x) 返回弧度x的余弦值。返回值范围为 [1,1]
    cosh(x) 双曲余弦(hyperbolic cosine)函数,计算x的双曲余弦值。
    cross(A,B) 返回两个三元向量的叉积(cross product)。注意,输入参数必须是三元向量!
    degrees(x) 输入参数为弧度值(radians),函数将其转换为角度值(degrees)
    determinant(m) 计算矩阵的行列式因子。
    dot(A,B) 返回AB的点积(dot product)。参数AB可以是标量,也可以是向量(输入参数方面,点积和叉积函数有很大不同)。
    exp(x) 计算 ex 的值,e=2.71828182845904523536
    exp2(x) 计算 2x 的值
    floor(x) 对输入参数向下取整。例如floor(float(1.3))返回的值为1.0;但是floor(float(-1.3))返回的值为-2.0。该函数与ceil(x)函数相对应。
    fmod(x,y) 返回x/y的余数。如果y0,结果不可预料。
    frac(x) 返回标量或矢量的小数
    frexp(x, out i) 将浮点数x分解为尾数和指数,即 x=m2i , 返回m,并将指数存入i中;如果x0,则尾数和指数都返回0
    isfinite(x) 判断标量或者向量中的每个数据是否是有限数,如果是返回true;否则返回false;
    isinf(x) 判断标量或者向量中的每个数据是否是无限,如果是返回true;否则返回false;
    isnan(x) 判断标量或者向量中的每个数据是否是非数据(not-a-number NaN),如果是返回true;否则返回false;
    ldexp(x, n) 计算 x2n 的值
    lerp(a, b, f) 计算 (1f)a+bf 或者 a+f(ba) 的值。即在下限a和上限b之间进行插值,f表示权值。注意,如果ab是向量,则权值f必须是标量或者等长的向量。
    lit(NdotL, NdotH, m) N表示法向量;
    L表示入射光向量;
    H表示半角向量;
    m表示高光系数。 
    函数计算环境光、散射光、镜面光的贡献,返回的4元向量。 
    X位表示环境光的贡献,总是1.0; 
    Y位代表散射光的贡献,如果  NL<0 ,则为0;否则为 NL  
    Z位代表镜面光的贡献,如果 NL<0  或者 NH<0 ,则位0;否则为 (NL)m ;
    W位始终位1.0
    log(x) 计算 ln(x) 的值,x必须大于0
    log2(x) 计算 log(x)2 的值,x必须大于0
    log10(x) 计算 log(x)10 的值,x必须大于0
    max(a, b) 比较两个标量或等长向量元素,返回最大值。
    min(a,b) 比较两个标量或等长向量元素,返回最小值。
    modf(x, out ip)x分解成整数和分数两部分,每部分都和x有着相同的符号,整数部分被保存在ip中,分数部分由函数返回
    mul(M, N) 矩阵M和矩阵N的积,计算方法如下
    mul(M,N)=M11M12M13M14M21M22M23M24M31M32M33M34M41M42M43M44N11N12N12N13N21N22N23N24N31N32N33N34N41N42N43N44
    mul(M, v) 矩阵M和列向量v的积,公式如下
    mul(M,v)=M11M12M13M14M21M22M23M24M31M32M33M34M41M42M43M44v1v2v3v4
    mul(v, M) 行向量v和矩阵M的积,公式如下
    mul(v,M)=[v1v2v3v4]M11M12M13M14M21M22M23M24M31M32M33M34M41M42M43M44
    noise(x) 根据它的参数类型,这个函数可以是一元、二元或三元噪音函数。返回的值在01之间,并且通常与给定的输入值一样
    pow(x, y) xy
    radians(x) 函数将角度值转换为弧度值
    round(x) 返回四舍五入值。
    rsqrt(x) x的平方根的倒数,x必须大于0
    saturate(x)x限制到[0,1]之间
    sign(x) 如果 x>0 则返回1;否则返回0
    sin(x) 输入参数为弧度,计算正弦值,返回值范围 为[-1,1]
    sincos(float x, out s, out c) 该函数是同时计算x的sin值和cos值,其中s=sin(x)c=cos(x)。该函数用于“同时需要计算sin值和cos值的情况”,比分别运算要快很多!
    sinh(x) 计算x的双曲正弦
    smoothstep(min, max, x)x位于minmax区间中。如果x=min,返回0;如果x=max,返回1;如果x在两者之间,按照下列公式返回数据:
    2(xminmaxmin)3+3(xminmaxmin)2
    step(a, x) 如果 x<a ,返回0;否则,返回1
    sqrt(x)x的平方根, x x必须大于0
    tan(x) 计算x正切值
    tanh(x) 计算x的双曲线切线
    transpose(M) 矩阵M的转置矩阵
    如果M是一个AxB矩阵,M的转置是一个BxA矩阵,它的第一列是M的第一行,第二列是M的第二行,第三列是M的第三行,等等


    (2)几何函数

    函数 功能描述
    distance(pt1, pt2) 两点之间的欧几里德距离(Euclidean distance)
    faceforward(N,I,Ng) 如果 NgI<0 ,返回N;否则返回-N
    length(v) 返回一个向量的模,即sqrt(dot(v,v))
    normalize(v) 返回v向量的单位向量
    reflect(I, N) 根据入射光纤方向I和表面法向量N计算反射向量,仅对三元向量有效
    refract(I,N,eta) 根据入射光线方向I,表面法向量N和折射相对系数eta,计算折射向量。如果对给定的eta,IN之间的角度太大,返回(0,0,0)。
    只对三元向量有效


    (3)纹理映射函数

    函数 功能描述
    tex1D(sampler1D tex, float s) 一维纹理查询
    tex1D(sampler1D tex, float s, float dsdx, float dsdy) 使用导数值(derivatives)查询一维纹理
    Tex1D(sampler1D tex, float2 sz) 一维纹理查询,并进行深度值比较
    Tex1D(sampler1D tex, float2 sz, float dsdx,float dsdy) 使用导数值(derivatives)查询一维纹理, 并进行深度值比较
    Tex1Dproj(sampler1D tex, float2 sq) 一维投影纹理查询
    Tex1Dproj(sampler1D tex, float3 szq) 一维投影纹理查询,并比较深度值
    Tex2D(sampler2D tex, float2 s) 二维纹理查询
    Tex2D(sampler2D tex, float2 s, float2 dsdx, float2 dsdy) 使用导数值(derivatives)查询二维纹理
    Tex2D(sampler2D tex, float3 sz) 二维纹理查询,并进行深度值比较
    Tex2D(sampler2D tex, float3 sz, float2 dsdx,float2 dsdy) 使用导数值(derivatives)查询二维纹理,并进行深度值比较
    Tex2Dproj(sampler2D tex, float3 sq) 二维投影纹理查询
    Tex2Dproj(sampler2D tex, float4 szq) 二维投影纹理查询,并进行深度值比较
    texRECT(samplerRECT tex, float2 s) 二维非投影矩形纹理查询(OpenGL独有)
    texRECT (samplerRECT tex, float3 sz, float2 dsdx,float2 dsdy) 二维非投影使用导数的矩形纹理查询(OpenGL独有)
    texRECT (samplerRECT tex, float3 sz) 二维非投影深度比较矩形纹理查询(OpenGL独有)
    texRECT (samplerRECT tex, float3 sz, float2 dsdx,float2 dsdy) 二维非投影深度比较并使用导数的矩形纹理查询(OpenGL独有)
    texRECT proj(samplerRECT tex, float3 sq) 二维投影矩形纹理查询(OpenGL独有)
    texRECT proj(samplerRECT tex, float3 szq) 二维投影矩形纹理深度比较查询(OpenGL独有)
    Tex3D(sampler3D tex, float s) 三维纹理查询
    Tex3D(sampler3D tex, float3 s, float3 dsdx, float3 dsdy) 结合导数值(derivatives)查询三维纹理
    Tex3Dproj(sampler3D tex, float4 szq) 查询三维投影纹理,并进行深度值比较
    texCUBE(samplerCUBE tex, float3 s) 查询立方体纹理
    texCUBE (samplerCUBE tex, float3 s, float3 dsdx, float3 dsdy) 结合导数值(derivatives)查询立方体纹理
    texCUBEproj (samplerCUBE tex, float4 sq) 查询投影立方体纹理


    在这个表中,每个函数第二个参数的名字指明了在执行纹理查询的时候,它的值是如果被使用的:

    • s表示这是一个一元、二元或三元纹理坐标。
    • z表示这是一个用来进行阴影贴图查找的深度比较值。
    • q表示这是一个透视值,在进行纹理查找之前,它被用来除以纹理坐标(s)。

    当你使用的纹理函数允许你指定一个深度比较值的时候,与之相关联的纹理单元必须被设置成深度比较纹理。否则,深度比较实际上不会被执行。 

    (4)偏导函数

    函数 功能描述
    ddx(a) 近似a关于屏幕空间x轴的偏导数
    ddy(a) 近似a关于屏幕空间y轴的偏导数


    (5)调试函数

    函数 功能描述
    void debug(float4 x) 如果在编译时设置了DEBUG,片段着 色程序中调用该函数可以将值x作为COLOR语义的最终输出;否则该函数什么也不做。
    更多相关内容
  • 总体CG方法是求解对称正定多右端大型线性方程组的非常有效的方法之一.利用矩阵Schur补的定义和性质给出总体CG方法残量的精确表达式.
  • 可以用来观看CG方法在计算均匀和非均匀特征值下的数值性态
  • 本文介绍了一种新的CG方法来求解非线性方程组。 该方法使用精确的线搜索实现了下降和全局收敛的条件。 与其他方法相比,在迭代次数和函数评估次数方面,数值结果都很好。
  • MATLAB实现CG算法

    2022-04-25 13:36:40
    MATLAB实现CG算法
  • 高等数值分析上机报告,高等数值分析课程设计报告,代码都是跑过的没有问题,matlab使用的是R2019A,有需要的小伙伴可以参考用。
  • cg法matlab代码概述 这是来自FEA软件手动填充的*.profiles 。 它还具有HTML格式(快速且正确地呈现来自先前提交的html表),(有点慢,Firefox出现问题,通常是当前提交),并且第一行和Feature列已修复,以便于表...
  • 支持LBFGS和CG等优化方法 CPU或GPU加速 Mapreduce并行化 梯度检查 易于配置 基线实验 #ACKNOWLEDGEMENTS包含的minFunc代码文件夹由Mark Schmidt()提供。 Quoc V. Le()提供了MATLAB Mapreduce。 #USAGE要运行代码...
  • 不等式约束-锥约束(线性,二阶锥和半定)的原始对偶内点法,锥约束的对数屏障方法 受限-以上各项的任意组合 开源的 根据2条款BSD许可发布 免费且可以使用开放源代码和封闭源代码的商业代码 多语言支持与C ++,...
  • 共轭梯度方法CG)MATLAB编程及其和Gauss_Seidel方法的一个比较 共轭梯度法是一种迭代法,不同于Jacobi,Gauss-Seidel和SOR方法。理论上最多n步找到真解,实际计算中考虑到舍入误差,一般迭代3n~5n步,每步的运算...

    共轭梯度方法(CG)MATLAB编程及其和Gauss_Seidel方法的一个比较

    共轭梯度法是一种迭代法,不同于Jacobi,Gauss-Seidel和SOR方法。理论上最多n步找到真解,实际计算中考虑到舍入误差,一般迭代3n~5n步,每步的运算量相当于举证乘向量的运算量。这种神奇的方法,对应稀疏矩阵特别有效。

    % 我们来造一个大型稀疏正定矩阵,比较共轭梯度法和Gauss_Seidei方法迭代速度上的差异。
    %为了方便,不妨构造一个三对角矩阵。
    clc
    clear
    profile on
    m =1000;%先来一个一千阶的方程组,便于CG方法和Gauss_Seidel方法的比较。
    e = ones (m,1);
    A = spdiags ([ e,4*e,e ],[-1,0,1],m,m);
    %full(A)
    true_x = rand(m,1);
    b = A * true_x;
    x0 = rand(m,1);
    eps = 1.0e-9;
    
    %测试一千万阶的方程组,CG方法在本机需要若干秒。此时,Gauss_Seidel方法卡死。
    [x_CG,~] = CG(A,b,x0,eps);
    error_CG = norm(x_CG - true_x,2);
    
    %来看看Gauss_Seidel方法的速度
    [x_G,~] = Gauss_Seidel(A,b,x0,eps);
    error_G = norm(x_G - true_x,2);
    
    % MATLAB对于特殊的一千万阶大规模方程组的求解也几乎是秒算,可见其矩阵运算功能的强大。
    %不愧矩阵实验室,毕竟是经过优化的。
    x_matlab = A\b;
    error_matlab = norm(x_matlab - true_x,2);
     profile viewer
     p = profile('info');
     profsave(p,'profile_results')

    Alt text
    可以看得出,CG方法在解这类问题上甩Gauss_Seidel不仅仅是用几条街来衡量的。这只是在一千阶的情况下。
    Alt text
    需要注意的是,Gauss_Seidel方法时间主要消耗在了用该方法矩阵A是否收敛的判断上。实际当中,这一步是多余的,因为我们可以通过若干步的迭代结果来判断是否收敛,一般不收敛,很快就会接待到无穷,很少有解震荡不收敛的情况产生。

    下面附上两个子函数:


    1、CG.m

    function [x,steps] = CG(A,b,x0,eps)
    r0 = b - A*x0;
    p0 = r0;
    if nargin == 3
        eps = 1.0e-6;
    end
    steps = 0;
    while 1
        if abs(p0) < eps
            break;
        end
        steps = steps + 1;
        a0 = r0'*r0/(p0'*A*p0);%多次用到可以存一步。
        x1 = x0 + a0*p0;
    
        r1 = r0 -a0*A*p0;
    
        b0 = r1'*r1/(r0'*r0);
        %这里的r'* r虽然后面可能还会用到,但是由于计算量不大,没有必要再设个新变量将
        %其存下了,内存上的开销划不来。
        p1 = r1 + b0*p0;
    
        %只是用到前后两层的向量,所以节省内存开销,计算完后面一层,可以往回覆盖掉没用
        %的变量。
    
        x0 = x1;
        r0 = r1;
        p0 = p1;
    
    end
    x = x0;
    end

    2、Gauss_Seidel.m

    function [x,out] = Gauss_Seidel(A,b,x0,eps,iterNum)
    %% 同Jacobi迭代,只是迭代函数f有所改变。
    n = length(b);
    if det(A) == 0
        error('No Unique Solution Or No Solution!');
    end
    for i = 1 : n
        j = 1;
        while A(i,i) == 0
            if A(j,i) ~= 0 && A(i,j) ~= 0
                tempV = A(j,:);
                A(j,:) = A(i,:);
                A(i,:) = tempV;
                tempS = b(j);
                b(j) = b(i);
                b(i) = tempS;
            else
                j = j + 1;
            end
            if j > n
                error('No Exchange Exist!');
            end
        end
    end
    D = diag(diag(A));
    invD = diag(1./diag(A));
    J = invD * (D - A);
    invDb = invD * b;
    H = eye(n) - tril(A)^-1*A;
    if max(abs(eig(H))) >= 1
        error('Gauss_Seidel Algorithm Cannot Convergence!')
    end
    
        function x = f(x)
            for l = 1:n
                temp_x = J(l,:)*x + invDb(l);
                x(l) = temp_x;
            end
        end
    x = x0;
    if nargin == 5
        for k = 1:iterNum-1
            x = f(x);
        end
        x_0 = x;
        x = f(x);
        out = norm(x-x_0,2);
    else
        if nargin == 3
            eps = 1.0e-6;
        end
        out = 0;
        while 1
            x_0 = x;
            x = f(x);
            out = out + 1 ;
            if norm( x - x_0 ,2 ) < eps
                break;
            end
        end
    end
    end

    对于实正定方程组而言,CG方法是收敛的。CG方法也有若干改进,比如用cholesky分解,这个后面会写。

    展开全文
  • 在本文中,CG方法的高级版本,即投影CG方法,用于ECT系统的图像重建。 将解空间投影到Krylov子空间中,并通过CG方法在低维特定子空间中解决反问题。 对气固两相流进行了静态和动态实验。 使用通过投影CG方法获得的...
  • cg法matlab代码GPUTUM:有限元求解器 GPUTUM FEM解算器是为解决FEM线性系统而编写的C ++ / CUDA库。 它旨在通过使用GPU硬件快速解决FEM系统。 该代码由美国盐湖城犹他大学科学计算与成像研究所的Zhisong Fu和T. ...
  • CG+BiCG+CGS.rar

    2020-09-04 11:01:17
    与传统的逐次过松弛(SOR)方法不同,CG算法没有自由参数可供选择。当矩阵A为非对称矩阵时,CG算法的一个直接推广就是所谓的双共轭梯度法(BiCG)。共轭梯度平方法(CGS)只需要乘以A,便能产生良好的近似解。
  • Hager和H.Zhang,算法851:CG_DESCENT,一种有保证下降的共轭梯度方法,数学软件中的ACM Transactions,32(2006),113-137。 [3] WW Hager和H. Zhang,《非线性共轭梯度法研究》,《太平洋最优化》,第2卷(2006年...
  • 行业资料-交通装置-CG发动机曲轴箱传动系统的装配方法.zip
  • cg法matlab代码代码_非凸面_PS 灵活的Matlba代码可实现强大的非凸光度学立体声,具有多种选择,包括: RGB或灰度图像 正交摄影机或透视摄影机 自阴影,投射阴影,高光等的鲁棒性可通过以下方法实现:简单的阈值化,...
  • 改良梯度下降法

    资料

    参考视频:
    【详细推导】【本视频还证明了收敛性】https://www.bilibili.com/video/BV16a4y1t76z?from=search&seid=35153938940534319
    【简单介绍】https://www.bilibili.com/video/BV1jE411u7bx?from=search&seid=35153938940534319

    梯度下降法

    本算法由梯度下降所引申:
    对于优化问题:
    在这里插入图片描述
    使用梯度下降:
    在这里插入图片描述
    注意,alpha也是可以算的:求一个a使得f(x0+aP0)最小,是求函数极值的,这时候是关于a的一个函数,所以对a求导求极小值,复合函数求导法则最后就会得到求导等于0,这时候的点是驻点,就是导数值为0的点,因为二阶导数黑塞矩阵正定,所以一定为极小值点。这时候就求出了在P0方向上的最小值点。
    图中()意味内积。

    共轭和预备知识

    共轭:A共轭的定义
    在这里插入图片描述
    引申:非零向量如果共轭,则一定线性无关。
    在这里插入图片描述
    此图:只要a_i彼此正交,则xi很容易计算。

    问题再次转化:要求x,可以用aplha表示,而alpha只和p有关
    因此,只要给出一组A共轭的向量组,就可以把x写出来。

    如何获得P?

    在这里插入图片描述
    怎么得到线性无关的向量组?
    施密特正交法:【不断地做A-投影(指内积是A内积)】
    在这里插入图片描述

    共轭梯度下降算法

    在这里插入图片描述
    残量是梯度的相反数】
    【图中的(AP0,R1)代表内积】
    一二步和最速下降一样,对新的残量。
    新的前进方向,让r1在P0上做投影,沿着共轭方向前进。

    优点

    • 可以证明,比梯度下降法收敛的快
    • 可以证明能得到精确解
      在这里插入图片描述
      图中span代表子空间,可以证明是同一个子空间,最后一个等号是科里洛夫子空间

    代码实现

    在这里插入图片描述
    在这里插入图片描述
    此图中,去掉min res的注释即可得到梯度下降法代码。

    例题

    在这里插入图片描述
    在这里插入图片描述
    第二步迭代:【经评论区老哥提醒,下图P1的正负号写反了,感谢】
    在这里插入图片描述

    展开全文
  • 基于Newton法改进的BFGS迭代算法与Newton-CG算法,侯麟,尚晓吉,本文主要研究了数值分析中数值优化与非线性方程组求解这两个重要问题。文中首先概述了数值优化与非线性方程组的关系,然后对BFGS�
  • 一种基于Cg语言在图形处理器GPU上实现加密的方法.pdf
  • CGProject:CG课程

    2021-07-06 20:06:12
    common.js 一些复用的工具方法。 style.css 程序界面样式定义。 pj1.html 多边形区域填充的程序入口。 pj1.js 多边形区域填充的程序控制。 polygon.js 画线、多边形区域填充算法实现。 pj2.html 3D立方体旋转的程序...
  • CG技术纹理

    2015-06-01 21:48:47
    数字CG技术制作纹理的经典教材,可以学到实用的动画特效纹理制作方法
  • 共轭梯度(CG)算法

    万次阅读 多人点赞 2017-11-16 14:27:38
    共轭梯度(CG方法计算数学与科学工程计算研究所 陆嵩简单介绍共轭梯度方法也是一种迭代方法,不同于Jacobi,Gauss-Seidel和SOR方法,理论上只要n步就能找到真解,实际计算中,考虑到舍入误差,一般迭代3n到5n步,...

    共轭梯度(CG)方法

    简单介绍

    共轭梯度方法也是一种迭代方法,不同于Jacobi,Gauss-Seidel和SOR方法,理论上只要n步就能找到真解,实际计算中,考虑到舍入误差,一般迭代3n到5n步,每步的运算量相当与矩阵乘向量的运算量,对稀疏矩阵特别有效。
    共轭梯度方法对于求解大型稀疏矩阵是很棒的方法,但是这个方法看起来总不是太靠谱。这个方法也不是越迭代精度越高,有时候可能迭代多了,反而出错,对迭代终止条件的选择,要求还是很高的。
    共轭梯度法收敛的快慢依赖于系数矩阵的谱分布情况,当特征值比较集中,系数矩阵的条件数很小,共轭梯度方法收敛得就快。“超线性收敛性”告诉我们,实际当中,我们往往需要更少的步数就能得到所需的精度的解。

    基本思想和原理核心

    考虑线性对称正定方程组: A x = b Ax = b Ax=b 可以转化为求解二次泛函 ϕ ( x ) = 1 2 x T A x − b T x 。 \phi(x) = \frac{1}{2} x^TAx-b^Tx。 ϕ(x)=21xTAxbTx
    的最小值问题,直接可以验证 ϕ ′ ( x ) = A x − b \phi&#x27;(x)=Ax-b ϕ(x)=Axb
    当然,选择其他的二次泛函也能满足这个条件导数等于这个的条件了,那是另外的一些方法了,这里不提。

    比如说这里取 A = d i a g ( 2 , 2 ) , b = [ 2 , 2 ] T A = diag(2,2),b=[2,2]^T A=diag(2,2),b=[2,2]T,那么就有 1 2 x T A x − b T x = x 1 2 + x 2 2 − 2 x 1 − 2 x 2 = ( x 1 − 1 ) 2 + ( x 2 − 1 ) 2 − 2 \frac{1}{2}x^TAx-bTx = x_1^2+x_2^2-2x_1-2x_2=(x_1-1)^2+(x_2-1)^2-2 21xTAxbTx=x12+x222x12x2=(x11)2+(x21)22
    这里可以看出它是一个尖朝下,经过源点的椭圆抛物面(当然,这个例子中椭圆是圆),且在 ( 1 , 1 ) (1,1) (1,1)点达到最小值-2。如下图所示(手绘的,比较丑,不过我已经很努力了):
    这里写图片描述

    在介绍共轭梯度法前,我们得先说一个这样一个二次泛函的性质,用大白话来说,就是这样:

    • 用一根xy平面上的直线穿过二次曲面在xy平面上的值形成的椭圆,那么这个直线上的二次泛函的值在直线被椭圆截出的线段的中点处取到最小。且这个最小值是唯一的,如图。
      这里写图片描述

    • 最小值的梯度方向垂直于上面提到的那根直线。

    • 最小值点位于位于这根直线的共轭超平面上(一维),不知道什么是共轭超平面,没有涉及到公式推导计算,这里可先不用管,只要知道它是由这根直线的法平面和 A A A以及真解有关的一个平面即可。

    这个其实很好理解,对于二次泛函形成的函数图像,我们沿着每个方向用一个平行于z轴的平面去切它,得到的是一个抛物面,放到坐标平面上,就是一根抛物线,抛物线的最小值点当然在中点处取到。注意,这里的负梯度方向,就是沿最速下降方向,它不一定能经过椭圆的圆心。

    以上是以最简单的二维的情况来说明这个问题,更高维的情况也如是。比如说三维的情况,我们可以想象有一个西瓜,它的内部每一点处有某一种性质,比如说甜度(或密度等等),假设它在正中心处甜度最高,最外面是瓜皮最不甜,依次往里甜度是呈抛物增长的。那么,我们拿一把刀,给它切一刀(不一定对半),会形成一个切面,那么这个切面是个椭圆,且在椭圆的正中间是最甜的(极值点),拿一根牙签,从这个最甜的点垂直切面往里面一戳,这个方向就是这一点负梯度方向。当然,牙签够长,你也不一定能正好西瓜正中心。因此才有了从最速下降法到共轭梯度法的改进。

    有了以上的思想,我们很容易就能推导出共轭梯度法的算法过程,不想摆公式,还是以切西瓜为例:

    • 首先我们拿一根足够长的牙签沿着西瓜任意位置将它穿透,那么前面的性质告诉我们没在西瓜中的牙签的中点就是这根牙签上的最值点。(这其实就是最速下降方法的步骤,最速下降法往往是沿着一个方向取到最小,步长以此来决定,接着换一个方向,同样一个过程……)

    • 接着找到最值点(牙签中点)的梯度方向,最速下降法就是以梯度方向接着找最小值点,但我们现在不这么做。梯度方向和原来的牙签的方向形成了一个面,我们试图在这个面里面找一个更好的方向。前面的性质告诉我们,这个切面的中点是这个面的最小值点,那么我就应该以牙签中点和这个切面的连线作为方向是最理想的。

    • 问题是这个切面的中点不好找,好在前面的性质告诉我们,这个切面的中点的梯度方向是垂直于这个切面的,把提到的这几个条件联立起来,其实很快就能找到牙签中点和切面中点的连线方向,以及我们所需要的步长。

    把以上的过程用公式摆开,就得到了共轭梯度法算法过程,如下算法框图所示:
    这里写图片描述

    当然,共轭梯度算法并不是很快,后来人们提出了很多方法去加快这个速度,比如说各种预优方法等等,都是后话,这里不提。

    程序代码和结果

    C代码

    #include<stdio.h>
    #include<stdlib.h>
    #include<math.h>
    #include<string.h>
    #define N 5
    #define epsilon  0.005
    
    void main()
    {
        void matrixTimesVec(double A[N][N], double b[N], double Ab[N]);
        double scalarProduct(double vec1[], double vec2[]);
        void vecPlus(double vec1[], double vec2[], double vec[]);
        void numPlusVec(double num, double vec0[], double vec[]);
        int i, j;
        static double A[N][N] = { 0 };
        static double b[N] = { 1, 1, 1, 1, 1 };
        static double x0[N] = { 1, 1, 1, 1, 1 };
        double x[N], r[N], p[N], w[N], alpha, rho00, rho0, rho1, beta;
        //生成一个大规模稀疏矩阵A,这里以三对角为例。
    
    
    
    
        for (i = 1; i < N - 1; i++)
        {
            A[i][i - 1] = 2;
            A[i][i] = 3;
            A[i][i + 1] = 1;
        }
        A[0][0] = 3; A[0][1] = 1;
        A[N - 1][N - 2] = 2; A[N - 1][N - 1] = 3;
    
    
        printf("\n要求解的示例方程组为:\n A ||| b ||| x0\n");
        for (i = 0; i < N; i++)
        {
            for (j = 0; j < N; j++)
            {
                printf("%f ", A[i][j]);
            }
            printf("||| %f||| %f\n", b[i], x0[i]);
        }
    
        //init
        memcpy(x, x0, N*sizeof(double));
        double Ax0[N];
        matrixTimesVec(A, x0, Ax0);
        numPlusVec(-1.0, Ax0, Ax0);
        vecPlus(b, Ax0, r);
        rho0 = scalarProduct(r, r);
        rho00 = rho0;
        memcpy(p, r, N*sizeof(double));
        do
        {
            matrixTimesVec(A, p, w);
            alpha = rho0 / (scalarProduct(p, w));
            double temp[N];
            numPlusVec(alpha, p, temp);
            double x_temp[N];
            vecPlus(x, temp, x_temp);
            memcpy(x, x_temp, N*sizeof(double));
    
            numPlusVec(-alpha, w, temp);
            double r_temp[N];
            vecPlus(r, temp, r_temp);
            memcpy(r, r_temp, N*sizeof(double));
    
            rho1 = scalarProduct(r, r);
    
            beta = rho1 / rho0;
            numPlusVec(beta, p, temp);
            vecPlus(r, temp, p);
            rho0 = rho1;
    
        } while (rho1 > epsilon);
    
        printf("\n\n");
        printf("方程组的解为:\n");
        for (i = 0; i < N; i++)
            printf("%f\n", x[i]);
        getchar();
    }
    
    void matrixTimesVec(double A[N][N], double b[N], double Ab[N])
    {
        int i, j;
    
    
        for (i = 0; i < N; i++)
        {
            Ab[i] = 0.0;
            for (j = 0; j < N; j++)
            {
                Ab[i] = Ab[i] + A[i][j] * b[j];
            }
        }
    }
    
    double scalarProduct(double vec1[], double vec2[])
    {
        double s = 0;
        int i;
        for (i = 0; i < N; i++)
        {
            s = s + vec1[i] * vec2[i];
        }
        return s;
    }
    void vecPlus(double vec1[], double vec2[], double vec[])
    {
        int i;
        for (i = 0; i < N; i++)
        {
            vec[i] = vec1[i] + vec2[i];
        }
    }
    void numPlusVec(double num, double vec0[], double vec[])
    {
        int i;
        for (i = 0; i < N; i++)
            vec[i] = num*vec0[i];
    
    }
    

    这里写图片描述

    MATLAB代码

    clc
    clear
    %% 共轭梯度算法并不是迭代的次数越多越精确,所以要合理地设置迭代终止的的条件。
    epsilon = 0.005;
    N = 5;
    A = zeros(N);
    for i=2:N-1
        A(i,i-1)=2;
        A(i,i) = 3;
        A(i,i+1) =1;
    end
    A(1,1)=3;
    A(1,2) = 1;
    A(N,N-1) = 2;
    A(N,N) = 3;
    b = ones(N,1);
    x0 = b;
    
    %% 初始化
    x = x0;
    r = b - A*x0;
    rho0 = conj(r')*r;
    rho00 = rho0;
    p=r;
    
    while(rho0>0.005)
        w= A*p;
        alpha = rho0/(conj(p')*w);
        x = x+alpha*p;
        r = r- alpha*w;
        rho1 = conj(r')*r;
        beta = rho1/rho0;
        p = r + beta*p;
        rho0 = rho1;
        x
    
    end
    A\b
    x
    
    展开全文
  • CG 空间的真实世界场景数据库设计 映射到虚拟 CG 空间的真实世界场景数据库设计 Takashi Tomii,工程学院,横滨国立大学,横滨,日本 240- 8501 Minako Kobayashi Ricoh Co. Ltd., Tokyo, Japan 107-8544 Hiroshi ...
  • USTC_CG:00106501

    2021-03-21 17:30:24
    中国科学技术大学《计算机图形学》课程...C ++面向对象编程思想和方法,了解基础的设计模式和架构思维 快速看懂及使用网上的代码,库及各种资料 从问题到抽象,到数学建模再到算法实现的方法 使用Unity3D引擎开发3D游戏
  • 资源分类:Python库 所属语言:Python 资源全名:cg-20.14.1.tar.gz 资源来源:官方 安装方法:https://lanzao.blog.csdn.net/article/details/101784059
  • 城市公共危机会导致重大人身伤亡与财产损失,科学的城市应急系统可有效应对和处理危机事件。信息技术的发展使得计算机能为城市应急救援提供科学决策支持。...CG与 GIS的城市应急系统的建设思路和方法
  • cg法matlab代码RSOpt(黎曼随机优化算法) 作者:,,, 最后页面更新时间:2019年5月31日 最新版本:1.0.3(有关更多信息,请参见发行说明) 审理 设f:M - > R上上一个M A光滑实值函数。 目标问题与M上的给定模型...
  • 模糊集matlab代码CG-nrIT2FCM 在这个项目中,我们专注于区间类型 2 模糊集生成和遥感影像分类。 用于 T1FCM、IT2FCM、s-IT2FCM 和 nr-IT2FCM 方法的 Matlab 软件。 作者:邢HH和于XC 电子邮件: 该软件实现了我论文...
  • 这种方法是基于Krylov子空间的方法中的一种 求解大型稀疏线性方程组
  • cg:talexander配置文件

    2021-02-15 21:28:39
    遵循此处记录的方法: : 设置 git init --bare $HOME /.cg alias cg= ' /usr/bin/git --git-dir=$HOME/.cg/ --work-tree=$HOME ' cg config --local status.showUntrackedFiles no 将文件添加到cg cg add ....

空空如也

空空如也

1 2 3 4 5 ... 20
收藏数 34,411
精华内容 13,764
关键字:

CG方法

友情链接: HW1.rar