精华内容
参与话题
问答
  • 高斯函数(Gaussian function)的详细分析

    万次阅读 多人点赞 2017-06-11 20:26:22
    ... 论文中遇到很重要的一个元素就是高斯核函数,但是必须要分析出高斯函数的各种潜在属性,本文首先参考相关材料给出高斯核函数的基础,然后使用matlab自动保存不同参数下的高斯核函数的变化gif动图,同时...

    原文地址:http://blog.csdn.net/jorg_zhao/article/details/52687448


    摘要

        论文中遇到很重要的一个元素就是高斯核函数,但是必须要分析出高斯函数的各种潜在属性,本文首先参考相关材料给出高斯核函数的基础,然后使用matlab自动保存不同参数下的高斯核函数的变化gif动图,同时分享出源代码,这样也便于后续的论文写作。


    高斯函数的基础

    2.1 一维高斯函数

    高斯函数,Gaussian Function, 也简称为Gaussian,一维形式如下:


    对于任意的实数a,b,c,是以著名数学家Carl Friedrich Gauss的名字命名的。高斯的一维图是特征对称“bell curve”形状,a是曲线尖峰的高度,b是尖峰中心的坐标,c称为标准方差,表征的是bell钟状的宽度。


    高斯函数广泛应用于统计学领域,用于表述正态分布,在信号处理领域,用于定义高斯滤波器,在图像处理领域,二维高斯核函数常用于高斯模糊Gaussian Blur,在数学领域,主要是用于解决热力方程和扩散方程,以及定义Weiertrass Transform。

    从上图可以看出,高斯函数是一个指数函数,其log函数是对数凹二次函数 whose logarithm a concave quadratic function。

    高斯函数的积分是误差函数error function,尽管如此,其在整个实线上的反常积分能够被精确的计算出来,使用如下的高斯积分


    同理可得


    当且仅当

    上式积分为1,在这种情况下,高斯是正态分布随机变量的概率密度函数,期望值μ=b,方差delta^2 = c^2,即



    2.2 二维高斯函数

        二维高斯函数,形如


    A是幅值,x。y。是中心点坐标,σσy是方差,图示如下,A = 1, xo = 0, yo = 0, σx = σy = 1



    2.3 高斯函数分析

    这一节使用matlab直观的查看高斯函数,在实际编程应用中,高斯函数中的参数有

    ksize 高斯函数的大小

    sigma 高斯函数的方差

    center 高斯函数尖峰中心点坐标

    bias 高斯函数尖峰中心点的偏移量,用于控制截断高斯函数

    为了方便直观的观察高斯函数参数改变而结果也不一样,下面的代码实现了参数的自动递增,并且将所有的结果图保存为gif图像,首先贴出完整代码:

    [plain] view plain copy
     print?
    1.  function mainfunc()  
    2. % 测试高斯函数,递增的方法实现高斯函数参数的改变对整个高斯函数的影响,  
    3. % 并自动保存为gif格式输出。  
    4. % created by zhao.buaa 2016.09.28  
    5.   
    6. %% 保存gif动画  
    7. item = 10;      % 迭代次数  
    8. dt = 1;             % 步长大小  
    9. ksize =20;      % 高斯大小  
    10. sigma = 2;      % 方差大小  
    11. % filename = ['ksize-' num2str(ksize) '--' num2str(ksize+dt*item) '-sigma-' num2str(sigma) '.gif']; %必须预先建立gif文件  
    12. filename = ['ksize-' num2str(ksize)  '-sigma-' num2str(sigma) '--' num2str(sigma+dt*item) '.gif']; %必须预先建立gif文件  
    13.   
    14. % main loop  
    15. for i = 1:item  
    16.     center  = round(ksize/2);          % 中心点  
    17.     bias       = ksize*10/10;              % 偏移中心点量  
    18.     ksigma = ksigma(ksize, sigma, center, bias);    % 构建核函数  
    19.     tname  = ['ksize-' num2str(ksize) '-sigma-' num2str(sigma) '-center-' num2str(center)];  
    20.     figure(i), mesh(ksigma), title(tname);  
    21.     %设置固定的x-y-z坐标范围,便于观察,axis([xmin xmax ymin ymax zmin zmax])  
    22.     axis([0 ksize 0 ksize 0 0.008]);  view([0, 90]);% 改变可视角度     
    23.     % ksize 递增  
    24. %     ksize = ksize + 10;  
    25.     % sigma 递增  
    26.     sigma = sigma + dt;       
    27.       
    28.     % 自动保存为gif图像  
    29.     frame = getframe(i);  
    30.     im = frame2im(frame);  
    31.     [I,map] = rgb2ind(im,256);  
    32.     if i==1  
    33.         imwrite(I,map,filename,'gif','Loopcount',inf, 'DelayTime',0.4);  
    34.     else  
    35.         imwrite(I,map,filename,'gif','WriteMode','append','DelayTime',0.4);  
    36.     end  
    37. end;  
    38.   
    39. close all;  
    40.   
    41.   
    42. %% 截断高斯核函数,截断的程度取决于参数bias  
    43. function ksigma = ksigma(ksize, sigma, center,bias)  
    44. %ksize = 80;    sigma = 15;  
    45. ksigma=fspecial('gaussian',ksize, sigma);   % 构建高斯函数  
    46. [m, n] =size(ksigma);  
    47. for i = 1:m  
    48.     for j = 1:n  
    49.         if(  (i<center-bias)||(i>center+bias)||(j<center-bias)||(j>center+bias)  )  
    50.             ksigma(i,j) = 0;  
    51.         end;  
    52.     end;  
    53. end;  

    结果图:

    固定ksize为20,sigma从1-9,固定center在高斯中间,并且bias偏移量为整个半径,即原始高斯函数。

    随着sigma的增大,整个高斯函数的尖峰逐渐减小,整体也变的更加平缓,则对图像的平滑效果越来越明显。

    保持参数不变,对上述高斯函数进行截断,即truncated gaussian function,bias的大小为ksize*3/10,则结果如下图:


    truncated gaussian function的作用主要是对超过一定区域的原始图像信息不再考虑,这就保证在更加合理的利用靠近高斯函数中心点的周围像素,同时还可以改变高斯函数的中心坐标,如下图:


    为了便于观察截断的效果,改变了可视角度。


    高斯核函数卷积

        论文中使用gaussian与feature map做卷积,目前的结果来看,要做到随着到边界的距离改变高斯函数的截断参数,因为图像的边缘如果使用原始高斯函数,就会在边界地方出现特别低的一圈,原因也很简单,高斯函数在与原始图像进行高斯卷积的时候,图像边缘外为0计算的,那么如何解决边缘问题呢?

    先看一段代码:

    [plain] view plain copy
     print?
    1. % 截断高斯核函数  
    2. ksize = 40; ksigma=fspecial('gaussian',  ksize, 6);  
    3. [m, n] =size(ksigma);  
    4. for i = 1:m  
    5.     for j = 1:n  
    6.         if( i<25 )  
    7.            ksigma(i,j) = 0;  
    8.         end;  
    9.     end;  
    10. end;  
    11. figure, mesh(ksigma);  

    在i,即row上对高斯核函数进行截断,bias为半径大小,则如图6

    并且对下图7进行卷积,

    首先卷积核为原始未截断高斯核函数,则结果如图8

    可以看出,在图像边缘处的卷积结果出现不想预见的结果,边缘处的值出现大幅度减少的情况,这是高斯核函数在边缘处将图像外的部分当成0计算的结果,因此,需要对高斯核函数进行针对性的截断处理,但是前提是要掌握bias的规律,下面就详细分析。

    如果使用图6的高斯核函数与图7做卷积操作,则如图9:

    可以看出,相比较于图8,与高斯核函数相对应的部分出现了变化,也就是说:

    [plain] view plain copy
     print?
    1. if( i<25 )  
    2.    ksigma(i,j) = 0;  
    3. end;  

    靠近边缘的时候,改变 i 或 j 的值,即可保证边缘处的平滑处理。但是这样改变高斯核函数,使用matlab不是很好解决这个问题,还是使用将待处理图像边缘向外部扩展bias的大小,与标准高斯核函数做卷积,再将超过原始图像大小的部分剪切掉,目前来看在使用matlab中imfilter函数做卷积运算最合适且最简单的处理方法了,先写在这里,此部分并不是论文中的核心部分,只是数值运算的技巧性编程方法。


    展开全文
  • Gauss

    2013-04-28 00:37:00
    1 //Guass_未测试版本 2 /* 3 如果可以肯定有唯一解,可以用化成上三角,时间有一点优势;但是如果存在自由变量 4 那么就不能这么写,而且化成上三角后,对应的行号和列号会发生错位; ... 9 void ga...
      1 //Guass_未测试版本
      2 /*
      3 如果可以肯定有唯一解,可以用化成上三角,时间有一点优势;但是如果存在自由变量
      4 那么就不能这么写,而且化成上三角后,对应的行号和列号会发生错位;
      5  
      6 */
      7 const double eps=1e-10;
      8 typedef double Matrix[N][N];
      9 void gauss(Matrix A,int n){ //n个方程组,n个自由变量;A[i][n]放常量y[i]; 
     10     int i,j,r,c,id;
     11     for (r=0;r<n;r++){
     12         id=r;
     13         for (i=r+1;i<n;i++) if (fabs(A[id][r])<fabs(A[i][r])) id=i;
     14         if (id!=r) for (j=r;j<=n;j++) swap(A[id][j],A[r][j]);
     15         
     16         for (i=r+1;i<n;i++){
     17             double tmp=A[i][r]/A[r][r];
     18             for (j=r;j<=n;j++){
     19                 A[i][j]-=tmp*A[r][j];
     20             } 
     21         }
     22     }
     23     //在A[i][n]放x[i]; 
     24     for (i=n-1;i>=0;i--){
     25         for (j=n-1;j>i;j--){
     26             A[i][n]-=A[i][j]*A[j][n];
     27         }
     28         A[i][n]/=A[i][i];
     29     }    
     30 }
     31 /*
     32 化成斜线方程,可以避免行号和列号发生错位,但时间会稍慢一点,
     33 而且当方程组有解时,唯一解:x[i]=A[i][n]/A[i][i]
     34 无穷多解,A[i][n]/A[i][i]表示在自由变量取0时的x[i]; 
     35 */ 
     36 int gauss_jordan(Matrix A,int n){
     37     int i,j,r,c,id;
     38     for (r=0;r<n;r++){
     39         id=r;
     40         for (i=r+1;i<n;i++) if (fabs(A[id][r])<fabs(A[i][r])) id=i;
     41         if (fabs(A[id][r])<eps) continue;
     42         
     43         for (i=0;i<n;i++){
     44             if (i==r) continue;
     45             double f=A[i][r]/A[r][r];
     46             for (j=r+1;j<=n;j++){
     47                 A[i][j]-=f*A[i][j];
     48             }
     49         }
     50     }
     51     int cnt=0;
     52     for (i=0;i<n;i++){
     53         if (fabs(A[i][i])<eps && fabs(A[i][n])>eps) return -1;//无解;
     54         if (fabs(A[i][i])<eps && fabs(A[i][n])<eps) cnt++;//自由变量; 
     55     }    
     56     if (cnt==0) return 0;//唯一解; 
     57     return cnt;//自由元的数量; 
     58 } 
     59 //提高精度版本;行列变换用它们的最小公倍数,所以只在最后求解的时候除法;
     60 //注意数据会不会越界;
     61 int gauss_jordan(Matrix A,int n){
     62     int i,j,r,c,id;
     63     for (r=0;r<n;r++){
     64         id=r;
     65         for (i=r+1;i<n;i++) if (abs(A[id][r])<abs(A[i][r])) id=i;
     66         if (abs(A[id][r])==0) continue;
     67         
     68         for (i=0;i<n;i++){
     69             if (i==r) continue;
     70             int g=__gcd(A[r][r],A[i][r]);
     71             int tp1=A[r][r]/g,tp2=A[i][r]/g;
     72             for (j=r;j<=n;j++){
     73                 A[i][j]=A[i][j]*tp1-A[r][j]*tp2;
     74             }
     75         }
     76     }
     77     int cnt=0;
     78     for (i=0;i<n;i++){
     79         if (abs(A[i][i])==0 && abs(A[i][n])>0) return -1;//无解;
     80         if (abs(A[i][i])==0 && abs(A[i][n])==0) cnt++;//自由变量; 
     81     }    
     82     if (cnt==0) return 0;//唯一解; 
     83     return cnt;//自由元的数量; 
     84 } 
     85 
     86 
     87 /*二进制的gauss*/
     88 int gauss(Matrix A,int n,int m){
     89     int i,j,r,c,id;
     90     for (r=0,c=0;r<n && c<m;r++,c++){
     91         
     92         for (i=r;i<n;i++) if (A[i][c]) { id=i;break; }
     93         if (A[id][c])==0){ r--;continue; }
     94         if (id!=r) for (j=c;j<=m;j++) swap(A[id][j],A[r][j]);
     95         
     96         for (i=r+1;i<n;i++)if (A[i][c]){
     97             for (j=c;j<=m;j++) A[i][j]^=A[r][j];
     98         }    
     99     }
    100     for (i=r;i<n;i++) if (abs(A[i][n])>0) return -1;//无解;
    101     if (m-r>0) return m-r;//自由元的个数;
    102      
    103 }

     hust 1651 **

    View Code
      1 /*
      2 题意:n*m的矩阵,要么是非负整数,要么是-1,-1表示要填的数,
      3 然后告诉你每一行,每一列的和,问你有无解,每行每列至多2个-1;
      4 
      5 很简单的想法就是列出n+m个方程组然后gauss,但gauss给我们的印象就是o(n^3;
      6 所以也许就把这个想法抛弃了,但实际是每个未知数最多出现在两个方程里,也就是
      7 在行里变换中,只需要变换一行后就可以直接break,时间复杂度就是O(n^2);
      8 
      9 */
     10 #include<iostream>
     11 #include<cmath>
     12 #include<vector>
     13 #include<algorithm>
     14 #include<cstdio>
     15 #include<cstdlib>
     16 #include<cstring>
     17 using namespace std;
     18 const int N=1000+10;
     19 typedef double Matrix[N][N];
     20 int mz[N][N],c[N][N];
     21 int a[N],b[N];
     22 int n,m,sz;
     23 const double eps=1e-10;
     24 Matrix A;
     25 int gauss(Matrix A,int n,int m){
     26     int i,j,r,c,id;
     27     for (r=0,c=0;r<n && c<m;r++,c++){
     28         id=r;
     29         for (i=r;i<n;i++) if (fabs(A[i][c])>fabs(A[id][c])) id=i;
     30         if (fabs(A[id][c])<eps) {  r--;continue; }
     31         if (id!=r) for (j=c;j<=m;j++) swap(A[id][j],A[r][j]);
     32 
     33         for (i=r+1;i<n;i++){
     34             if (A[i][c]==0) continue;
     35             double f=A[i][c]/A[r][c];
     36             for (j=c;j<=m;j++){
     37                 A[i][j]-=f*A[r][j];
     38             }
     39             break;
     40         }
     41     }
     42     for (i=r;i<n;i++) if (abs(A[i][m])>0) return -1;
     43     if (m-r>0) return 2;
     44     return 1;
     45 }
     46 void work(){
     47     memset(A,0,sizeof(A));
     48     for (int i=0;i<n;i++){
     49         int sum=0;
     50         for (int j=0;j<m;j++){
     51             if (c[i][j]!=-1) A[i][c[i][j]]=1;
     52             sum+=mz[i][j];
     53         }
     54         A[i][sz]=a[i]-sum;
     55     }
     56     for (int j=0;j<m;j++){
     57         int sum=0;
     58         for (int i=0;i<n;i++){
     59             if (c[i][j]!=-1) A[j+n][c[i][j]]=1;
     60             sum+=mz[i][j];
     61         }
     62         A[n+j][sz]=b[j]-sum;
     63     }
     64   /*  cout<<"++++ "<<endl;
     65     for (int i=0;i<n+m;i++){
     66         for (int j=0;j<=sz;j++){
     67             cout<<A[i][j]<<" ";
     68         }cout<<endl;
     69     }
     70     cout<<"+++++ "<<endl;
     71     */int t=gauss(A,n+m,sz);
     72     if (t==-1) printf("No solution\n");
     73     else if (t==1) printf("Unique\n");
     74     else printf("More than one\n");
     75 }
     76 int main(){
     77     int T,cas=0;
     78     scanf("%d",&T);
     79     while (T--){
     80         scanf("%d%d",&n,&m);
     81         sz=0;
     82         memset(c,-1,sizeof(c));
     83         for (int i=0;i<n;i++){
     84             for (int j=0;j<m;j++){
     85                 scanf("%d",&mz[i][j]);
     86                 if (mz[i][j]==-1){
     87                     mz[i][j]=0;
     88                     c[i][j]=sz++;
     89                 }
     90             }
     91         }
     92         for (int i=0;i<n;i++) scanf("%d",&a[i]);
     93         for (int j=0;j<m;j++) scanf("%d",&b[j]);
     94         printf("Case #%d: ",++cas);
     95         work();
     96 
     97 
     98     }
     99 
    100     return 0;
    101 }
    102 /*
    103 5
    104 2 2
    105 1 1
    106 -1 -1
    107 2 4
    108 3 3
    109 2 2
    110 1 1
    111 -1 -1
    112 2 5
    113 3 3
    114 2 2
    115 1 1
    116 2 2
    117 2 4
    118 3 3
    119 4 4
    120 1 2 3 4
    121 1 -1 -1 2
    122 3 -1 -1 4
    123 5 6 7 8
    124 10 6 14 26
    125 10 13 15 18
    126 2 2
    127 1 -1
    128 2 -1
    129 -1 1
    130 3 -3
    131 
    132 */

     

    转载于:https://www.cnblogs.com/Rlemon/archive/2013/04/28/3048344.html

    展开全文
  • Matlab计算Jacobi、Gauss-Seidel迭代例题Jacobi迭代Gauss-Seidel迭代 例题 方程组{5x1+2x2+x3=−12−x1+4x2+2x3=202x1−3x2+10x3=3\begin{cases} 5x_1 + 2x_2 + x_3 = -12\\ -x_1 + 4x_2 + 2x_3 = 20\\ 2x_1 - 3x_2 ...

    使用Matlab实现:Jacobi、Gauss-Seidel迭代

    例题

    方程组{5x1+2x2+x3=12x1+4x2+2x3=202x13x2+10x3=3\begin{cases} 5x_1 + 2x_2 + x_3 = -12\\ -x_1 + 4x_2 + 2x_3 = 20\\ 2x_1 - 3x_2 + 10x_3 = 3\\ \end{cases} 求解,当maxxi(k+1)xi(k)105max|x_i^{(k + 1)} - x_i^{(k)}| \leq 10^{-5} 时候迭代终止。

    以下解答过程,上标表示迭代次数,下标表示序号。

    Jacobi迭代

    定义变量:
    D=diag(a11,a22,...,ann),L=[0a210...ai1...ai,i10...an1...an,i1...an,n10],U=[0a12...a1,i...a1,n...0ai1,i...ai1,n...0an1,n...0]D = diag(a_{11}, a_{22}, ..., a_{nn}),\\ L = \left[\begin{array}{cccccc} 0\\ -a_{21} &amp; 0\\ ...\\ -a_{i1} &amp; ... &amp; -a_{i,i-1} &amp; 0\\ ...\\ -a_{n1} &amp; ... &amp; -a_{n,i-1} &amp; ... &amp; -a_{n,n-1} &amp; 0\\ \end{array}\right],\\ U = \left[\begin{array}{cccccc} 0 &amp;-a_{12} &amp; ... &amp; -a_{1,i} &amp; ... &amp; -a_{1,n}\\ ...\\ &amp; &amp; 0 &amp; -a_{i-1,i} &amp; ... &amp; -a_{i-1,n}\\ ...\\ &amp; &amp; &amp; &amp; 0 &amp; -a_{n-1,n}\\ ...\\ &amp; &amp; &amp; &amp; &amp; 0\\ \end{array}\right]

    其矩阵迭代形式为:
    x(k+1)=BJx(k)+fJBJ=D1(L+U),fJ=D1bx^{(k+1)} = B_J \cdot x^{(k)} + f_J\\ B_J = D^{-1} \cdot (L + U), \quad f_J = D^{-1} \cdot b

    写出分量形式:

    {x1(k+1)=15(122x2(k)x3(k))x2(k+1)=14(20+x1(k)2x3(k))x3(k+1)=110(32x1(k)+3x2(k))\begin{cases} x_1^{(k + 1)} = \frac15 (-12 - 2 x_2^{(k)} - x_3^{(k)} )\\ x_2^{(k + 1)} = \frac14 (20 + x_1^{(k)} - 2x_3^{(k)} )\\ x_3^{(k + 1)} = \frac1{10} (3 - 2 x_1^{(k)} + 3x_2^{(k)} )\\ \end{cases}

    写成矩阵形式:

    [x1(k+1)x2(k+1)x3(k+1)]=[00.40.20.2500.50.20.30][x1(k)x2(k)x3(k)]+[2.450.3]\left[\begin{array}{c} x_1^{(k + 1)}\\ x_2^{(k + 1)}\\ x_3^{(k + 1)}\\ \end{array}\right] = \left[\begin{array}{cccc} 0 &amp; -0.4 &amp; -0.2\\ 0.25 &amp; 0 &amp; -0.5\\ -0.2 &amp; 0.3 &amp; 0\\ \end{array}\right] \cdot \left[\begin{array}{c} x_1^{(k)}\\ x_2^{(k)}\\ x_3^{(k)}\\ \end{array}\right] + \left[\begin{array}{c} -2.4\\ 5\\ 0.3\\ \end{array}\right]

    取初始向量:x0=(0,0,0)Tx^0 = (0, 0, 0)^T,依次按照上式进行迭代。使用Matlab进行编程求解。

    a=[0,-0.4,-0.2;0.25,0,-0.5;-0.2,0.3,0];
    b = [-2.4;5;0.3];
    x = [0;0;0];
    xx = a * x + b;
    i = 0;
    while norm(x - xx, inf) >= 1e-5
    	x = xx;
    	xx = a * x + b;
    	i = i +1;
    end
    

    以上代码,最终 x=xi,xx=x(i+1)x = x^{i}, xx = x^{(i + 1)} ,最终迭代次数位 i+1i + 1 次,如果你需要看到更长的小数位置,可以使用以下Matlab代码,表示使用15位浮点或定点数。

    format long g
    

    运行结果为:
    在这里插入图片描述

    即精确解为 x=(4,3,2)Tx = (-4,3,2)^T

    Gauss-Seidel迭代

    其矩阵迭代形式为:
    x(k+1)=BGx(k)+fGBG=(DL)1U,fG=(DL)1bx^{(k+1)} = B_G \cdot x^{(k)} + f_G\\ B_G = (D - L) ^{-1} \cdot U, \quad f_G = (D - L) ^{-1} \cdot b

    使用Matlab编程求解:

    d = [5,0,0;0,4,0;0,0,10];
    l = [0,0,0;1,0,0;-2,3,0];
    u = [0,-2,-1;0,0,-2;0,0,0];
    b = [-12;20;3];
    t = inv(d - l);
    bg = t * u;
    fg = t * b;
    x = [0;0;0];
    xx = [-2.4;4.4;2.1];
    i = 0;
    while norm(x - xx, inf) >= 1e-5
    	x = xx;
    	xx = bg * x + fg;
    	i = i +1;
    end
    

    运行结果为:
    在这里插入图片描述

    同样求得精确解为 x=(4,3,2)Tx = (-4,3,2)^T

    展开全文
  • 数值优化之高斯-牛顿法(Gauss-Newton)

    万次阅读 多人点赞 2018-04-12 19:51:48
    我的机器学习教程「美团」算法工程师带你入门机器学习 已经开始更新了,欢迎大家订阅~ 任何关于算法、编程、AI行业知识或博客内容的问题,可以随时扫码关注公众号「图灵的猫」,加入”学习小组“,沙雕博主在线答疑...

    我的机器学习教程「美团」算法工程师带你入门机器学习   已经开始更新了,欢迎大家订阅~

    任何关于算法、编程、AI行业知识或博客内容的问题,可以随时扫码关注公众号「图灵的猫」,加入”学习小组“,沙雕博主在线答疑~此外,公众号内还有更多AI、算法、编程和大数据知识分享,以及免费的SSR节点和学习资料。其他平台(知乎/B站)也是同名「图灵的猫」,不要迷路哦~

     

     

     

     

     

    一、基本概念定义

     

    1.非线性方程定义及最优化方法简述

       指因变量与自变量之间的关系不是线性的关系,比如平方关系、对数关系、指数关系、三角函数关系等等。对于此类方程,求解n元实函数f在整个n维向量空间Rn上的最优值点往往很难得到精确解,经常需要求近似解问题。

       求解该最优化问题的方法大多是逐次一维搜索的迭代算法,基本思想是在一个近似点处选定一个有利于搜索方向,沿这个方向进行一维搜索,得到新的近似点。如此反复迭代,知道满足预定的精度要求为止。根据搜索方向的取法不同,这类迭代算法可分为两类:

    解析法:需要用目标函数的到函数,

    梯度法:又称最速下降法,是早期的解析法,收敛速度较慢

    牛顿法:收敛速度快,但不稳定,计算也较困难。高斯牛顿法基于其改进,但目标作用不同

    共轭梯度法:收敛较快,效果好

    变尺度法:效率较高,常用DFP法(Davidon Fletcher Powell)

     

     

    直接法:不涉及导数,只用到函数值。有交替方向法(又称坐标轮换法)、模式搜索法、旋转方向法、鲍威尔共轭方向法和单纯形加速法等。

     

    2.非线性最小二乘问题

       非线性最小二乘问题来自于非线性回归,即通过观察自变量和因变量数据,求非线性目标函数的系数参数,使得函数模型与观测量尽量相似。

       高斯牛顿法解决非线性最小二乘问题的最基本方法,并且它只能处理二次函数。(使用时必须将目标函数转化为二次的)

       Unlike Newton'smethod, the Gauss–Newton algorithm can only be used to minimize a sum ofsquared function values

     

    3.基本数学表达

    a.梯度gradient,由多元函数的各个偏导数组成的向量

    以二元函数为例,其梯度为:

     

    b.黑森矩阵Hessian matrix,由多元函数的二阶偏导数组成的方阵,描述函数的局部曲率,以二元函数为例,

     

    c.雅可比矩阵 Jacobian matrix,是多元函数一阶偏导数以一定方式排列成的矩阵,体现了一个可微方程与给出点的最优线性逼近。以二元函数为例,

     

    如果扩展多维的话F: Rn-> Rm,则雅可比矩阵是一个m行n列的矩阵:

     

     

    雅可比矩阵作用,如果P是Rn中的一点,F在P点可微分,那么在这一点的导数由JF(P)给出,在此情况下,由F(P)描述的线性算子即接近点P的F的最优线性逼近:

     

    d.残差 residual,表示实际观测值与估计值(拟合值)之间的差

     

     

    二、牛顿法

    牛顿法的基本思想是采用多项式函数来逼近给定的函数值,然后求出极小点的估计值,重复操作,直到达到一定精度为止。

    1.考虑如下一维无约束的极小化问题:

     

    因此,一维牛顿法的计算步骤如下:

     

     

     

    需要注意的是,牛顿法在求极值的时候,如果初始点选取不好,则可能不收敛于极小点

     

    2.下面给出多维无约束极值的情形:

    若非线性目标函数f(x)具有二阶连续偏导,在x(k)为其极小点的某一近似,在这一点取f(x)的二阶泰勒展开,即:

     

      如果f(x)是二次函数,则其黑森矩阵H为常数,式(1)是精确的(等于号),在这种情况下,从任意一点处罚,用式(2)只要一步可求出f(x)的极小点(假设黑森矩阵正定,所有特征值大于0)

      如果f(x)不是二次函数,式(1)仅是一个近似表达式,此时,按式(2)求得的极小点,只是f(x)的近似极小点。在这种情况下,常按照下面选取搜索方向:

    牛顿法收敛的速度很快,当f(x)的二阶导数及其黑森矩阵的逆矩阵便于计算时,这一方法非常有效。【但通常黑森矩阵很不好求】

     

    3.下面给出一个实际计算例子。

     

    例:试用牛顿法求的极小值

     

    解:

     

     

    【f(x)是二次函数,H矩阵为常数,只要任意点出发,只要一步即可求出极小点】

     

    三、牛顿高斯法

    1.      gauss-newton是如何由上述派生的

    有时候为了拟合数据,比如根据重投影误差求相机位姿(R,T为方程系数),常常将求解模型转化为非线性最小二乘问题。高斯牛顿法正是用于解决非线性最小二乘问题,达到数据拟合、参数估计和函数估计的目的。

    假设我们研究如下形式的非线性最小二乘问题:

     

    这两个位置间残差(重投影误差):

     

    如果有大量观测点(多维),我们可以通过选择合理的T使得残差的平方和最小求得两个相机之间的位姿。机器视觉这块暂时不扩展,接着说怎么求非线性最小二乘问题。

    若用牛顿法求式3,则牛顿迭代公式为:

     

    看到这里大家都明白高斯牛顿和牛顿法的差异了吧,就在这迭代项上。经典高斯牛顿算法迭代步长λ为1.

    那回过头来,高斯牛顿法里为啥要舍弃黑森矩阵的二阶偏导数呢?主要问题是因为牛顿法中Hessian矩阵中的二阶信息项通常难以计算或者花费的工作量很大,而利用整个H的割线近似也不可取,因为在计算梯度时已经得到J(x),这样H中的一阶信息项JTJ几乎是现成的。鉴于此,为了简化计算,获得有效算法,我们可用一阶导数信息逼近二阶信息项。注意这么干的前提是,残差r接近于零或者接近线性函数从而接近与零时,二阶信息项才可以忽略。通常称为“小残量问题”,否则高斯牛顿法不收敛。

     

     

    3.  举例

    接下来的代码里并没有保证算法收敛的机制,在例子2的自嗨中可以看到劣势。关于自变量维数,代码可以支持多元,但两个例子都是一维的,比如例子1中只有年份t,其实可以增加其他因素的,不必在意。

     

    例子1,根据美国1815年至1885年数据,估计人口模型中的参数A和B。如下表所示,已知年份和人口总量,及人口模型方程,求方程中的参数。

     

     

     

    // A simple demo of Gauss-Newton algorithm on a user defined function  
      
    #include <cstdio>  
    #include <vector>  
    #include <opencv2/core/core.hpp>  
      
    using namespace std;  
    using namespace cv;  
      
    const double DERIV_STEP = 1e-5;  
    const int MAX_ITER = 100;  
      
      
    void GaussNewton(double(*Func)(const Mat &input, const Mat ¶ms), // function pointer  
                     const Mat &inputs, const Mat &outputs, Mat ¶ms);  
      
    double Deriv(double(*Func)(const Mat &input, const Mat ¶ms), // function pointer  
                 const Mat &input, const Mat ¶ms, int n);  
      
    // The user defines their function here  
    double Func(const Mat &input, const Mat ¶ms);  
      
    int main()  
    {  
        // For this demo we're going to try and fit to the function  
        // F = A*exp(t*B)  
        // There are 2 parameters: A B  
        int num_params = 2;  
      
        // Generate random data using these parameters  
        int total_data = 8;  
      
        Mat inputs(total_data, 1, CV_64F);  
        Mat outputs(total_data, 1, CV_64F);  
      
        //load observation data  
        for(int i=0; i < total_data; i++) {  
            inputs.at<double>(i,0) = i+1;  //load year  
        }  
        //load America population  
        outputs.at<double>(0,0)= 8.3;  
        outputs.at<double>(1,0)= 11.0;  
        outputs.at<double>(2,0)= 14.7;  
        outputs.at<double>(3,0)= 19.7;  
        outputs.at<double>(4,0)= 26.7;  
        outputs.at<double>(5,0)= 35.2;  
        outputs.at<double>(6,0)= 44.4;  
        outputs.at<double>(7,0)= 55.9;  
      
        // Guess the parameters, it should be close to the true value, else it can fail for very sensitive functions!  
        Mat params(num_params, 1, CV_64F);  
      
        //init guess  
        params.at<double>(0,0) = 6;  
        params.at<double>(1,0) = 0.3;  
      
        GaussNewton(Func, inputs, outputs, params);  
      
        printf("Parameters from GaussNewton: %f %f\n", params.at<double>(0,0), params.at<double>(1,0));  
      
        return 0;  
    }  
      
    double Func(const Mat &input, const Mat ¶ms)  
    {  
        // Assumes input is a single row matrix  
        // Assumes params is a column matrix  
      
        double A = params.at<double>(0,0);  
        double B = params.at<double>(1,0);  
      
        double x = input.at<double>(0,0);  
      
        return A*exp(x*B);  
    }  
      
    //calc the n-th params' partial derivation , the params are our  final target  
    double Deriv(double(*Func)(const Mat &input, const Mat ¶ms), const Mat &input, const Mat ¶ms, int n)  
    {  
        // Assumes input is a single row matrix  
      
        // Returns the derivative of the nth parameter  
        Mat params1 = params.clone();  
        Mat params2 = params.clone();  
      
        // Use central difference  to get derivative  
        params1.at<double>(n,0) -= DERIV_STEP;  
        params2.at<double>(n,0) += DERIV_STEP;  
      
        double p1 = Func(input, params1);  
        double p2 = Func(input, params2);  
      
        double d = (p2 - p1) / (2*DERIV_STEP);  
      
        return d;  
    }  
      
    void GaussNewton(double(*Func)(const Mat &input, const Mat ¶ms),  
                     const Mat &inputs, const Mat &outputs, Mat ¶ms)  
    {  
        int m = inputs.rows;  
        int n = inputs.cols;  
        int num_params = params.rows;  
      
        Mat r(m, 1, CV_64F); // residual matrix  
        Mat Jf(m, num_params, CV_64F); // Jacobian of Func()  
        Mat input(1, n, CV_64F); // single row input  
      
        double last_mse = 0;  
      
        for(int i=0; i < MAX_ITER; i++) {  
            double mse = 0;  
      
            for(int j=0; j < m; j++) {  
                for(int k=0; k < n; k++) {//copy Independent variable vector, the year  
                    input.at<double>(0,k) = inputs.at<double>(j,k);  
                }  
      
                r.at<double>(j,0) = outputs.at<double>(j,0) - Func(input, params);//diff between estimate and observation population  
      
                mse += r.at<double>(j,0)*r.at<double>(j,0);  
      
                for(int k=0; k < num_params; k++) {  
                    Jf.at<double>(j,k) = Deriv(Func, input, params, k);  
                }  
            }  
      
            mse /= m;  
      
            // The difference in mse is very small, so quit  
            if(fabs(mse - last_mse) < 1e-8) {  
                break;  
            }  
      
            Mat delta = ((Jf.t()*Jf)).inv() * Jf.t()*r;  
            params += delta;  
      
            //printf("%d: mse=%f\n", i, mse);  
            printf("%d %f\n", i, mse);  
      
            last_mse = mse;  
        }  
    }  
    
     
     


    运行结果:

     

     

    A=7.0,B=0.26  (初始值,A=6,B=0.3),100次迭代到第4次就收敛了。

    若初始值A=1,B=1,则要迭代14次收敛。

     

    下图为根据上面得到的A、B系数,利用matlab拟合的人口模型曲线

     

    例子2:我想要拟合如下模型,

     

    由于缺乏观测量,就自导自演,假设4个参数已知A=5,B=1,C=10,D=2,构造100个随机数作为x的观测值,计算相应的函数观测值。然后,利用这些观测值,反推4个参数。

     

     

     

    // A simple demo of Gauss-Newton algorithm on a user defined function  
      
    #include <cstdio>  
    #include <vector>  
    #include <opencv2/core/core.hpp>  
      
    using namespace std;  
    using namespace cv;  
      
    const double DERIV_STEP = 1e-5;  
    const int MAX_ITER = 100;  
      
      
    void GaussNewton(double(*Func)(const Mat &input, const Mat ¶ms), // function pointer  
                     const Mat &inputs, const Mat &outputs, Mat ¶ms);  
      
    double Deriv(double(*Func)(const Mat &input, const Mat ¶ms), // function pointer  
                 const Mat &input, const Mat ¶ms, int n);  
      
    // The user defines their function here  
    double Func(const Mat &input, const Mat ¶ms);  
      
    int main()  
    {  
        // For this demo we're going to try and fit to the function  
        // F = A*sin(Bx) + C*cos(Dx)  
        // There are 4 parameters: A, B, C, D  
        int num_params = 4;  
      
        // Generate random data using these parameters  
        int total_data = 100;  
      
        double A = 5;  
        double B = 1;  
        double C = 10;  
        double D = 2;  
      
        Mat inputs(total_data, 1, CV_64F);  
        Mat outputs(total_data, 1, CV_64F);  
      
        for(int i=0; i < total_data; i++) {  
            double x = -10.0 + 20.0* rand() / (1.0 + RAND_MAX); // random between [-10 and 10]  
            double y = A*sin(B*x) + C*cos(D*x);  
      
            // Add some noise  
           // y += -1.0 + 2.0*rand() / (1.0 + RAND_MAX);  
      
            inputs.at<double>(i,0) = x;  
            outputs.at<double>(i,0) = y;  
        }  
      
        // Guess the parameters, it should be close to the true value, else it can fail for very sensitive functions!  
        Mat params(num_params, 1, CV_64F);  
      
        params.at<double>(0,0) = 1;  
        params.at<double>(1,0) = 1;  
        params.at<double>(2,0) = 8; // changing to 1 will cause it not to find the solution, too far away  
        params.at<double>(3,0) = 1;  
      
        GaussNewton(Func, inputs, outputs, params);  
      
        printf("True parameters: %f %f %f %f\n", A, B, C, D);  
        printf("Parameters from GaussNewton: %f %f %f %f\n", params.at<double>(0,0), params.at<double>(1,0),  
                                                            params.at<double>(2,0), params.at<double>(3,0));  
      
        return 0;  
    }  
      
    double Func(const Mat &input, const Mat ¶ms)  
    {  
        // Assumes input is a single row matrix  
        // Assumes params is a column matrix  
      
        double A = params.at<double>(0,0);  
        double B = params.at<double>(1,0);  
        double C = params.at<double>(2,0);  
        double D = params.at<double>(3,0);  
      
        double x = input.at<double>(0,0);  
      
        return A*sin(B*x) + C*cos(D*x);  
    }  
      
    double Deriv(double(*Func)(const Mat &input, const Mat ¶ms), const Mat &input, const Mat ¶ms, int n)  
    {  
        // Assumes input is a single row matrix  
      
        // Returns the derivative of the nth parameter  
        Mat params1 = params.clone();  
        Mat params2 = params.clone();  
      
        // Use central difference  to get derivative  
        params1.at<double>(n,0) -= DERIV_STEP;  
        params2.at<double>(n,0) += DERIV_STEP;  
      
        double p1 = Func(input, params1);  
        double p2 = Func(input, params2);  
      
        double d = (p2 - p1) / (2*DERIV_STEP);  
      
        return d;  
    }  
      
    void GaussNewton(double(*Func)(const Mat &input, const Mat ¶ms),  
                     const Mat &inputs, const Mat &outputs, Mat ¶ms)  
    {  
        int m = inputs.rows;  
        int n = inputs.cols;  
        int num_params = params.rows;  
      
        Mat r(m, 1, CV_64F); // residual matrix  
        Mat Jf(m, num_params, CV_64F); // Jacobian of Func()  
        Mat input(1, n, CV_64F); // single row input  
      
        double last_mse = 0;  
      
        for(int i=0; i < MAX_ITER; i++) {  
            double mse = 0;  
      
            for(int j=0; j < m; j++) {  
                for(int k=0; k < n; k++) {  
                    input.at<double>(0,k) = inputs.at<double>(j,k);  
                }  
      
                r.at<double>(j,0) = outputs.at<double>(j,0) - Func(input, params);  
      
                mse += r.at<double>(j,0)*r.at<double>(j,0);  
      
                for(int k=0; k < num_params; k++) {  
                    Jf.at<double>(j,k) = Deriv(Func, input, params, k);  
                }  
            }  
      
            mse /= m;  
      
            // The difference in mse is very small, so quit  
            if(fabs(mse - last_mse) < 1e-8) {  
                break;  
            }  
      
            Mat delta = ((Jf.t()*Jf)).inv() * Jf.t()*r;  
            params += delta;  
      
            //printf("%d: mse=%f\n", i, mse);  
            printf("%f\n",mse);  
      
            last_mse = mse;  
        }  
    }  
     
     


    运行结果,得到的参数并不够理想,50次后收敛了

     

    下图中,每次迭代残差并没有持续减少,有反复

     

    4.优缺点分析

    优点:

    • 对于零残量问题,即r=0,有局部二阶收敛速度
    • 对于小残量问题,即r较小或接近线性,有快的局部收敛速度
    • 对于线性最小二乘问题,一步达到极小点

     

    缺点:

    • 对于不是很严重的大残量问题,有较慢的局部收敛速度
    • 对于残量很大的问题或r的非线性程度很大的问题,不收敛
    • 不一定总体收敛
    • 如果J不满秩,则方法无定义

     

    对于它的缺点,我们通过增加线性搜索策略,保证目标函数每一步下降,对于几乎所有非线性最小二乘问题,它都具有局部收敛性及总体收敛,即所谓的阻尼高斯牛顿法。

     

    其中,ak为一维搜索因子

    展开全文
  • 文章目录初始GaussDBGaussDB的版本GaussDB版本的区别OLTP和OLAP比较GaussDB T介绍GaussDB A 介绍MPP架构介绍架构组件介绍 初始GaussDB 名字的由来:GaussDB是华为数据库产品品牌名,致敬数据加高斯(GaussGaussDB...
  • GaussView5.08.rar是高斯软件的代表作品,这款是正式的应用版而且是免费的!!!
  • 通过GaussView建立模型,以便于在Gaussian中计算过渡态能量和结构、键和反应能量、分子轨道、原子电荷和电势、振动频率、红外和拉曼光谱、核磁性质、极化率和超极化率、热力学性质、反应路径,计算可以对体系的基态...
  • 题目:编写Gauss求积法计算积分的程序(Gauss点数取1,2,3,4,5即可)并用于计算积分 KaTeX parse error: No such environment: equation at position 8: \begin{̲e̲q̲u̲a̲t̲i̲o̲n̲}̲I=\int_{0}^{1} … 部分...
  • 高斯函数(Gaussian function)的详细分析

    万次阅读 多人点赞 2018-09-03 14:49:49
    摘要  论文中遇到很重要的一个元素就是高斯核函数,但是必须要分析出高斯函数的各种潜在属性,本文首先参考相关材料给出高斯核函数的基础,然后使用matlab自动保存不同参数下的高斯核函数的变化gif动图,同时分享...
  • Jacobi迭代法与Gauss-Seidel迭代法

    万次阅读 多人点赞 2016-01-25 19:29:58
    今天刚好有朋友和我讨论泊松图像融合算法,我说我过去文章里给出的是最原始、最直观的实现算法。对于理解泊松融合的原理比较有帮助,但是效率可能并不理想。印象中,泊松融合是有一个以矩阵为基础的快速算法的。...
  • GaussDB数据库Linux安装

    万次阅读 2019-08-14 15:25:20
    一、GaussDB: 全球首款AI-Native数据库 二、FusionStorage 8.0:业界性能第一的分布式存储 三、下载安装包 四、安装步骤 五、华为GaussDB数据库知识扩展 一、GaussDB: 全球首款AI-Native数据库 作为全球首款...
  • GaussView5

    2017-10-26 20:38:03
    Gview是一个专门设计与高斯配套使用的软件,其主要用途有两个:构建高斯的输入文件和以图的形式显示高斯计算的结果。除了可以自己构建输入文件外,Gview还可读入CHEM3D,...现在比较常用的有GaussView3和GaussView5.
  • Gauss-Newton算法学习

    万次阅读 多人点赞 2016-06-08 20:11:08
    Gauss-Newton算法是解决非线性最优问题的常见算法之一,最近研读DPPTAM开源项目代码,又碰到了,索性深入看下。本次讲解内容如下:基本数学名词识记牛顿法推导、算法步骤、计算实例高斯牛顿法推导(如何从牛顿法派生)...
  • Linux安装GaussDB数据库图文

    万次阅读 2019-08-23 14:20:00
    目录 一、说明 二、环境搭建 三、下载链接 四、官方安装文档目录(详情见上方链接) ...第一种版本适合大部分企业和开发人员,X86服务器+CentOS7/Red Hat 第二种版本适合大部分国产化的企业,鲲鹏服务器+麒麟操作...
  • Gauss-Seidel迭代法

    2012-12-15 21:46:54
    Gauss-Seidel迭代法的基本思想 由Jacobi迭代法中,每一次的迭代只用到前一次的迭代值,若每一次迭代充分利用当前最新的迭代值,即在计算第 个分量 时,用最新分量 , 代替旧分量 , ,就得到所谓解方程组的Gauss-...
  • GaussDB系列数据库简介

    2020-08-07 10:45:56
    GaussDB目前支持TD、Natezza、Oracle、MySQL、DB2、sybase、PG的离线数据迁移, 支持Oracle的全量+增强的在线迁移。 GaussDB版本的区别: GaussDB T(OLTP): 前身是GauussDB 100, 主打OLTP在线事务处理。 用于存储...
  • gaussview 5.08

    2015-11-21 08:25:41
    gaussview是gaussian的可视化窗口,在gaussview中可构建分子结构并选择合适的算法直接提交gaussian进行计算。
  • 1.创建数据库表和索引(数据库操作工具使用的是华为Data Studio) 设置数据库主键列: 创建系列,用于生成唯一主键ID: ...•下载GaussDB数据库JDBC驱动jar包: JDBC包名: com.huawei.gaus...
  • 华为GaussDB

    千次阅读 2019-09-23 14:02:44
    华为于2019年5月15日在北京召开了数据库及存储产品发布会,发布了人工智能原生(AI-Native)数据库GaussDB和分布式存储FusionStorage 8.0共两款产品。发布会上华为常务董事汪涛提出:“华为此次面向全球发布AI-...
  • Jacobi迭代、Gauss_Seidel迭代和最佳因子SOR迭代的比较 我们知道无论是有限差分方法还是有限元方法,最后都离不开解大规模方程组。除提到的这三种方法外,常用的方法还有共轭梯度法(CG)和多重网格方法以及他们的...

空空如也

1 2 3 4 5 ... 20
收藏数 21,172
精华内容 8,468
关键字:

gauss