最小二乘法–高斯牛顿迭代法

2023-11-18

最小二乘法–高斯牛顿迭代法

本文将详解最小二乘法的非线性拟合,高斯牛顿迭代法。

1.原理

高斯—牛顿迭代法的基本思想是使用泰勒级数展开式去近似地代替非线性回归模型,然后通过多次迭代,多次修正回归系数,使回归系数不断逼近非线性回归模型的最佳回归系数,最后使原模型的残差平方和达到最小。

①已知m个点:

这里写图片描述

②函数原型:

这里写图片描述

其中:(m>=n)

这里写图片描述

③目的是找到最优解β,使得残差平方和最小:

这里写图片描述

残差:

这里写图片描述

④要求最小值,即S的对β偏导数等于0:

这里写图片描述

⑤在非线性系统中,这里写图片描述是变量和参数的函数,没有close解。因此,我们给定一个初始值,用迭代法逼近解:

这里写图片描述

其中k是迭代次数,这里写图片描述是迭代矢量。

⑥而每次迭代函数是线性的,在这里写图片描述处用泰勒级数展开:

这里写图片描述

其中:J是已知的矩阵,为了方便迭代,令这里写图片描述

⑦此时残差表示为:

这里写图片描述

这里写图片描述

⑧带入公式④有:

这里写图片描述

化解得:

这里写图片描述

⑨写成矩阵形式:

这里写图片描述

⑩所以最终迭代公式为:

这里写图片描述

其中,Jf是函数f=(x,β)对β的雅可比矩阵。

2.代码

Java代码实现,解维基百科的例子:

https://en.wikipedia.org/wiki/Gauss%E2%80%93Newton_algorithm

①已知数据:

这里写图片描述

②函数模型:

这里写图片描述

③残差公式:

这里写图片描述

④对β求偏导数:

这里写图片描述

⑤代码如下:

<code class="hljs cs has-numbering"><span class="hljs-keyword">public</span> <span class="hljs-keyword">class</span> GaussNewton {
    <span class="hljs-keyword">double</span>[] xData = <span class="hljs-keyword">new</span> <span class="hljs-keyword">double</span>[]{<span class="hljs-number">0.038</span>, <span class="hljs-number">0.194</span>, <span class="hljs-number">0.425</span>, <span class="hljs-number">0.626</span>, <span class="hljs-number">1.253</span>, <span class="hljs-number">2.500</span>, <span class="hljs-number">3.740</span>};
    <span class="hljs-keyword">double</span>[] yData = <span class="hljs-keyword">new</span> <span class="hljs-keyword">double</span>[]{<span class="hljs-number">0.050</span>, <span class="hljs-number">0.127</span>, <span class="hljs-number">0.094</span>, <span class="hljs-number">0.2122</span>, <span class="hljs-number">0.2729</span>, <span class="hljs-number">0.2665</span>, <span class="hljs-number">0.3317</span>};

    <span class="hljs-keyword">double</span>[][] bMatrix = <span class="hljs-keyword">new</span> <span class="hljs-keyword">double</span>[<span class="hljs-number">2</span>][<span class="hljs-number">1</span>];<span class="hljs-comment">//系数β矩阵</span>
    <span class="hljs-keyword">int</span> m = xData.length;
    <span class="hljs-keyword">int</span> n = bMatrix.length;
    <span class="hljs-keyword">int</span> iterations = <span class="hljs-number">7</span>;<span class="hljs-comment">//迭代次数</span>

    <span class="hljs-comment">//迭代公式求解,即1中公式⑩</span>
    <span class="hljs-keyword">private</span> <span class="hljs-keyword">void</span> <span class="hljs-title">magic</span>(){
        <span class="hljs-comment">//β1,β2迭代初值</span>
        bMatrix[<span class="hljs-number">0</span>][<span class="hljs-number">0</span>] = <span class="hljs-number">0.9</span>;
        bMatrix[<span class="hljs-number">1</span>][<span class="hljs-number">0</span>] = <span class="hljs-number">0.2</span>;

        <span class="hljs-comment">//求J矩阵</span>
        <span class="hljs-keyword">for</span>(<span class="hljs-keyword">int</span> k = <span class="hljs-number">0</span>; k < iterations; k++) { 
            <span class="hljs-keyword">double</span>[][] J = <span class="hljs-keyword">new</span> <span class="hljs-keyword">double</span>[m][n];
            <span class="hljs-keyword">double</span>[][] JT = <span class="hljs-keyword">new</span> <span class="hljs-keyword">double</span>[n][m];
            <span class="hljs-keyword">for</span>(<span class="hljs-keyword">int</span> i = <span class="hljs-number">0</span>; i < m; i++){
                <span class="hljs-keyword">for</span>(<span class="hljs-keyword">int</span> j = <span class="hljs-number">0</span>; j < n; j++) {
                    J[i][j] = secondDerivative(xData[i], bMatrix[<span class="hljs-number">0</span>][<span class="hljs-number">0</span>], bMatrix[<span class="hljs-number">1</span>][<span class="hljs-number">0</span>], j);
                }
            }

            JT = MatrixMath.invert(J);<span class="hljs-comment">//求转置矩阵JT</span>
            <span class="hljs-keyword">double</span>[][] invertedPart = MatrixMath.mult(JT, J);<span class="hljs-comment">//矩阵JT与J相乘</span>

            <span class="hljs-comment">//矩阵invertedPart行列式的值:|JT*J|</span>
            <span class="hljs-keyword">double</span> result = MatrixMath.mathDeterminantCalculation(invertedPart);

            <span class="hljs-comment">//求矩阵invertedPart的逆矩阵:(JT*J)^-1</span>
            <span class="hljs-keyword">double</span>[][] reversedPart = MatrixMath.getInverseMatrix(invertedPart, result);

            <span class="hljs-comment">//求方差r(β)矩阵: ri = yi - f(xi, b)</span>
            <span class="hljs-keyword">double</span>[][] residuals = <span class="hljs-keyword">new</span> <span class="hljs-keyword">double</span>[m][<span class="hljs-number">1</span>];
            <span class="hljs-keyword">for</span>(<span class="hljs-keyword">int</span> i = <span class="hljs-number">0</span>; i < m; i++) {
                residuals[i][<span class="hljs-number">0</span>] = yData[i] - (bMatrix[<span class="hljs-number">0</span>][<span class="hljs-number">0</span>] * xData[i]) / (bMatrix[<span class="hljs-number">1</span>][<span class="hljs-number">0</span>] + xData[i]);
            }

            <span class="hljs-comment">//求矩阵积reversedPart*JT*residuals: (JT*J)^-1*JT*r</span>
            <span class="hljs-keyword">double</span>[][] product = MatrixMath.mult(MatrixMath.mult(reversedPart, JT), residuals);

            <span class="hljs-comment">//迭代公式, 即公式⑩</span>
            <span class="hljs-keyword">double</span>[][] newB = MatrixMath.plus(bMatrix, product);
            bMatrix = newB;        
        }        
        <span class="hljs-comment">//显示系数值</span>
        System.<span class="hljs-keyword">out</span>.println(<span class="hljs-string">"b1: "</span> + bMatrix[<span class="hljs-number">0</span>][<span class="hljs-number">0</span>] + <span class="hljs-string">"\nb2: "</span> + bMatrix[<span class="hljs-number">1</span>][<span class="hljs-number">0</span>]);        
    }

    <span class="hljs-comment">//2中公式④</span>
    <span class="hljs-keyword">private</span> <span class="hljs-keyword">static</span> <span class="hljs-keyword">double</span> <span class="hljs-title">secondDerivative</span>(<span class="hljs-keyword">double</span> x, <span class="hljs-keyword">double</span> b1, <span class="hljs-keyword">double</span> b2, <span class="hljs-keyword">int</span> index){
        <span class="hljs-keyword">switch</span>(index) {
            <span class="hljs-keyword">case</span> <span class="hljs-number">0</span>: <span class="hljs-keyword">return</span> x / (b2 + x);<span class="hljs-comment">//对系数bi求导</span>
            <span class="hljs-keyword">case</span> <span class="hljs-number">1</span>: <span class="hljs-keyword">return</span> -<span class="hljs-number">1</span> * (b1 * x) / Math.pow((b2+x), <span class="hljs-number">2</span>);<span class="hljs-comment">//对系数b2求导</span>
        }
        <span class="hljs-keyword">return</span> <span class="hljs-number">0</span>;
    }

    <span class="hljs-keyword">public</span> <span class="hljs-keyword">static</span> <span class="hljs-keyword">void</span> <span class="hljs-title">main</span>(String[] args) {
        GaussNewton app = <span class="hljs-keyword">new</span> GaussNewton();
        app.magic();        
    }
}</code><ul class="pre-numbering" style="opacity: 0.454786;"><li>1</li><li>2</li><li>3</li><li>4</li><li>5</li><li>6</li><li>7</li><li>8</li><li>9</li><li>10</li><li>11</li><li>12</li><li>13</li><li>14</li><li>15</li><li>16</li><li>17</li><li>18</li><li>19</li><li>20</li><li>21</li><li>22</li><li>23</li><li>24</li><li>25</li><li>26</li><li>27</li><li>28</li><li>29</li><li>30</li><li>31</li><li>32</li><li>33</li><li>34</li><li>35</li><li>36</li><li>37</li><li>38</li><li>39</li><li>40</li><li>41</li><li>42</li><li>43</li><li>44</li><li>45</li><li>46</li><li>47</li><li>48</li><li>49</li><li>50</li><li>51</li><li>52</li><li>53</li><li>54</li><li>55</li><li>56</li><li>57</li><li>58</li><li>59</li><li>60</li><li>61</li><li>62</li><li>63</li><li>64</li><li>65</li></ul><div class="save_code tracking-ad" style="display: none;" data-mod="popu_249"><a target=_blank target="_blank"><img src="http://static.blog.csdn.net/images/save_snippets.png" alt="" /></a></div><ul class="pre-numbering"><li>1</li><li>2</li><li>3</li><li>4</li><li>5</li><li>6</li><li>7</li><li>8</li><li>9</li><li>10</li><li>11</li><li>12</li><li>13</li><li>14</li><li>15</li><li>16</li><li>17</li><li>18</li><li>19</li><li>20</li><li>21</li><li>22</li><li>23</li><li>24</li><li>25</li><li>26</li><li>27</li><li>28</li><li>29</li><li>30</li><li>31</li><li>32</li><li>33</li><li>34</li><li>35</li><li>36</li><li>37</li><li>38</li><li>39</li><li>40</li><li>41</li><li>42</li><li>43</li><li>44</li><li>45</li><li>46</li><li>47</li><li>48</li><li>49</li><li>50</li><li>51</li><li>52</li><li>53</li><li>54</li><li>55</li><li>56</li><li>57</li><li>58</li><li>59</li><li>60</li><li>61</li><li>62</li><li>63</li><li>64</li><li>65</li></ul>

运行,输出得到:

b1: 0.3618366954234483
b2: 0.5562654497238557

这里写图片描述

⑥其中用到的矩阵运算代码如下:

<code class="hljs java has-numbering"><span class="hljs-keyword">public</span> <span class="hljs-class"><span class="hljs-keyword">class</span> <span class="hljs-title">MatrixMath</span> {</span>

    <span class="hljs-javadoc">/**
     * 矩阵基本运算:加、减、乘、除、转置
     */</span>
        <span class="hljs-keyword">public</span> <span class="hljs-keyword">final</span> <span class="hljs-keyword">static</span> <span class="hljs-keyword">int</span> OPERATION_ADD = <span class="hljs-number">1</span>;  
        <span class="hljs-keyword">public</span> <span class="hljs-keyword">final</span> <span class="hljs-keyword">static</span> <span class="hljs-keyword">int</span> OPERATION_SUB = <span class="hljs-number">2</span>;  
        <span class="hljs-keyword">public</span> <span class="hljs-keyword">final</span> <span class="hljs-keyword">static</span> <span class="hljs-keyword">int</span> OPERATION_MUL = <span class="hljs-number">4</span>;        

        <span class="hljs-javadoc">/**
         * 矩阵加法
         *<span class="hljs-javadoctag"> @param</span> a 加数
         *<span class="hljs-javadoctag"> @param</span> b 被加数
         *<span class="hljs-javadoctag"> @return</span> 和
         */</span>
        <span class="hljs-keyword">public</span> <span class="hljs-keyword">static</span> <span class="hljs-keyword">double</span>[][] <span class="hljs-title">plus</span>(<span class="hljs-keyword">double</span>[][] a, <span class="hljs-keyword">double</span>[][] b){
            <span class="hljs-keyword">if</span>(legalOperation(a, b, OPERATION_ADD)) { 
                <span class="hljs-keyword">for</span>(<span class="hljs-keyword">int</span> i=<span class="hljs-number">0</span>; i<a.length; i++) {         
                    <span class="hljs-keyword">for</span>(<span class="hljs-keyword">int</span> j=<span class="hljs-number">0</span>; j<a[<span class="hljs-number">0</span>].length; j++) {  
                        a[i][j] = a[i][j] + b[i][j];  
                    }
                }   
            }
            <span class="hljs-keyword">return</span> a;
        }

        <span class="hljs-javadoc">/**
         * 矩阵减法
         *<span class="hljs-javadoctag"> @param</span> a 减数
         *<span class="hljs-javadoctag"> @param</span> b 被减数
         *<span class="hljs-javadoctag"> @return</span> 差
         */</span>
        <span class="hljs-keyword">public</span> <span class="hljs-keyword">static</span> <span class="hljs-keyword">double</span>[][] <span class="hljs-title">substract</span>(<span class="hljs-keyword">double</span>[][] a, <span class="hljs-keyword">double</span>[][] b){
            <span class="hljs-keyword">if</span>(legalOperation(a, b, OPERATION_SUB)) { 
                <span class="hljs-keyword">for</span>(<span class="hljs-keyword">int</span> i=<span class="hljs-number">0</span>; i<a.length; i++) {  
                    <span class="hljs-keyword">for</span>(<span class="hljs-keyword">int</span> j=<span class="hljs-number">0</span>; j<a[<span class="hljs-number">0</span>].length; j++) {  
                        a[i][j] = a[i][j] - b[i][j];  
                    }
                }       
            }
            <span class="hljs-keyword">return</span> a;
        }       

        <span class="hljs-javadoc">/**
         * 判断矩阵行列是否符合运算
         *<span class="hljs-javadoctag"> @param</span> a 进行运算的矩阵
         *<span class="hljs-javadoctag"> @param</span> b 进行运算的矩阵
         *<span class="hljs-javadoctag"> @param</span> type 运算类型
         *<span class="hljs-javadoctag"> @return</span> 符合/不符合运算
         */</span>
        <span class="hljs-keyword">private</span> <span class="hljs-keyword">static</span> <span class="hljs-keyword">boolean</span> <span class="hljs-title">legalOperation</span>(<span class="hljs-keyword">double</span>[][] a, <span class="hljs-keyword">double</span>[][] b, <span class="hljs-keyword">int</span> type) {  
            <span class="hljs-keyword">boolean</span> legal = <span class="hljs-keyword">true</span>;  
            <span class="hljs-keyword">if</span>(type == OPERATION_ADD || type == OPERATION_SUB)  
            {  
                <span class="hljs-keyword">if</span>(a.length != b.length || a[<span class="hljs-number">0</span>].length != b[<span class="hljs-number">0</span>].length) {  
                    legal = <span class="hljs-keyword">false</span>;  
                }  
            }   
            <span class="hljs-keyword">else</span> <span class="hljs-keyword">if</span>(type == OPERATION_MUL)  
            {  
                <span class="hljs-keyword">if</span>(a[<span class="hljs-number">0</span>].length != b.length) {  
                    legal = <span class="hljs-keyword">false</span>;  
                }  
            }  
            <span class="hljs-keyword">return</span> legal;  
        } 

        <span class="hljs-javadoc">/**
         * 矩阵乘法
         *<span class="hljs-javadoctag"> @param</span> a 乘数
         *<span class="hljs-javadoctag"> @param</span> b 被乘数
         *<span class="hljs-javadoctag"> @return</span> 积
         */</span>
        <span class="hljs-keyword">public</span> <span class="hljs-keyword">static</span> <span class="hljs-keyword">double</span>[][] <span class="hljs-title">mult</span>(<span class="hljs-keyword">double</span>[][] a, <span class="hljs-keyword">double</span>[][] b){
            <span class="hljs-keyword">if</span>(legalOperation(a, b, OPERATION_MUL)) {
                <span class="hljs-keyword">double</span>[][] result = <span class="hljs-keyword">new</span> <span class="hljs-keyword">double</span>[a.length][b[<span class="hljs-number">0</span>].length];
                <span class="hljs-keyword">for</span>(<span class="hljs-keyword">int</span> i=<span class="hljs-number">0</span>; i< a.length; i++) {  
                    <span class="hljs-keyword">for</span>(<span class="hljs-keyword">int</span> j=<span class="hljs-number">0</span>; j< b[<span class="hljs-number">0</span>].length; j++) {  
                        result[i][j] = calculateSingleResult(a, b, i, j);  
                    }
                }
                <span class="hljs-keyword">return</span> result;
            }
            <span class="hljs-keyword">else</span>
            {
                <span class="hljs-keyword">return</span> <span class="hljs-keyword">null</span>;
            }       
        }

        <span class="hljs-javadoc">/**
         * 矩阵乘法
         *<span class="hljs-javadoctag"> @param</span> a 乘数
         *<span class="hljs-javadoctag"> @param</span> b 被乘数
         *<span class="hljs-javadoctag"> @return</span> 积
         */</span>
        <span class="hljs-keyword">public</span> <span class="hljs-keyword">static</span> <span class="hljs-keyword">double</span>[][] <span class="hljs-title">mult</span>(<span class="hljs-keyword">double</span>[][] a, <span class="hljs-keyword">int</span> b) {  
            <span class="hljs-keyword">for</span>(<span class="hljs-keyword">int</span> i=<span class="hljs-number">0</span>; i<a.length; i++) {  
                <span class="hljs-keyword">for</span>(<span class="hljs-keyword">int</span> j=<span class="hljs-number">0</span>; j<a[<span class="hljs-number">0</span>].length; j++) {  
                    a[i][j] = a[i][j] * b;  
                }  
            }  
            <span class="hljs-keyword">return</span> a;  
        }

        <span class="hljs-javadoc">/**
         * 乘数矩阵的行元素与被乘数矩阵列元素积的和
         *<span class="hljs-javadoctag"> @param</span> a 乘数矩阵
         *<span class="hljs-javadoctag"> @param</span> b 被乘数矩阵
         *<span class="hljs-javadoctag"> @param</span> row 行
         *<span class="hljs-javadoctag"> @param</span> column 列
         *<span class="hljs-javadoctag"> @return</span> 值
         */</span>
        <span class="hljs-keyword">private</span> <span class="hljs-keyword">static</span> <span class="hljs-keyword">double</span> <span class="hljs-title">calculateSingleResult</span>(<span class="hljs-keyword">double</span>[][] a, <span class="hljs-keyword">double</span>[][] b, <span class="hljs-keyword">int</span> row, <span class="hljs-keyword">int</span> column) {  
            <span class="hljs-keyword">double</span> result = <span class="hljs-number">0.0</span>;  
            <span class="hljs-keyword">for</span>(<span class="hljs-keyword">int</span> k = <span class="hljs-number">0</span>; k< a[<span class="hljs-number">0</span>].length; k++) {  
                result += a[row][k] * b[k][column];  
            }  
            <span class="hljs-keyword">return</span> result;  
        }  

        <span class="hljs-javadoc">/**
         * 矩阵的转置
         *<span class="hljs-javadoctag"> @param</span> a 要转置的矩阵
         *<span class="hljs-javadoctag"> @return</span> 转置后的矩阵
         */</span>
        <span class="hljs-keyword">public</span> <span class="hljs-keyword">static</span> <span class="hljs-keyword">double</span>[][] <span class="hljs-title">invert</span>(<span class="hljs-keyword">double</span>[][] a){
            <span class="hljs-keyword">double</span>[][] result = <span class="hljs-keyword">new</span> <span class="hljs-keyword">double</span>[a[<span class="hljs-number">0</span>].length][a.length];
            <span class="hljs-keyword">for</span>(<span class="hljs-keyword">int</span> i=<span class="hljs-number">0</span>;i<a.length;i++){
                <span class="hljs-keyword">for</span>(<span class="hljs-keyword">int</span> j=<span class="hljs-number">0</span>;j<a[<span class="hljs-number">0</span>].length;j++){  
                      result[j][i]=a[i][j];  
                  }  
              }  
            <span class="hljs-keyword">return</span> result;
        }   

    <span class="hljs-javadoc">/** 
     * 求可逆矩阵(使用代数余子式的形式) 
     */</span>   
        <span class="hljs-javadoc">/** 
         * 求传入的矩阵的逆矩阵 
         *<span class="hljs-javadoctag"> @param</span> value 需要转换的矩阵 
         *<span class="hljs-javadoctag"> @return</span> 逆矩阵 
         */</span>  
        <span class="hljs-keyword">public</span> <span class="hljs-keyword">static</span> <span class="hljs-keyword">double</span>[][] <span class="hljs-title">getInverseMatrix</span>(<span class="hljs-keyword">double</span>[][] value,<span class="hljs-keyword">double</span> result){  
            <span class="hljs-keyword">double</span>[][] transferMatrix = <span class="hljs-keyword">new</span> <span class="hljs-keyword">double</span>[value.length][value[<span class="hljs-number">0</span>].length];  
            <span class="hljs-comment">//计算代数余子式,并赋值给|A|  </span>
            <span class="hljs-keyword">for</span> (<span class="hljs-keyword">int</span> i = <span class="hljs-number">0</span>; i < value.length; i++) {  
                <span class="hljs-keyword">for</span> (<span class="hljs-keyword">int</span> j = <span class="hljs-number">0</span>; j < value[i].length; j++) {  
                    transferMatrix[j][i] =  mathDeterminantCalculation(getNewMatrix(i, j, value));  
                    <span class="hljs-keyword">if</span> ((i+j)%<span class="hljs-number">2</span>!=<span class="hljs-number">0</span>) {  
                        transferMatrix[j][i] = -transferMatrix[j][i];  
                    }  
                    transferMatrix[j][i] /= result;   
                }  
            }  
            <span class="hljs-keyword">return</span> transferMatrix;  
        }  

        <span class="hljs-javadoc">/*** 
         * 求行列式的值
         *<span class="hljs-javadoctag"> @param</span> value 需要算的行列式 
         *<span class="hljs-javadoctag"> @return</span> 计算的结果 
         */</span>  
        <span class="hljs-keyword">public</span> <span class="hljs-keyword">static</span> <span class="hljs-keyword">double</span> <span class="hljs-title">mathDeterminantCalculation</span>(<span class="hljs-keyword">double</span>[][] value){  
            <span class="hljs-keyword">if</span> (value.length == <span class="hljs-number">1</span>) {  
                <span class="hljs-comment">//当行列式为1阶的时候就直接返回本身  </span>
                <span class="hljs-keyword">return</span> value[<span class="hljs-number">0</span>][<span class="hljs-number">0</span>];  
            }

            <span class="hljs-keyword">if</span> (value.length == <span class="hljs-number">2</span>) {  
                <span class="hljs-comment">//如果行列式为二阶的时候直接进行计算  </span>
                <span class="hljs-keyword">return</span> value[<span class="hljs-number">0</span>][<span class="hljs-number">0</span>]*value[<span class="hljs-number">1</span>][<span class="hljs-number">1</span>]-value[<span class="hljs-number">0</span>][<span class="hljs-number">1</span>]*value[<span class="hljs-number">1</span>][<span class="hljs-number">0</span>];  
            }  

            <span class="hljs-comment">//当行列式的阶数大于2时  </span>
            <span class="hljs-keyword">double</span> result = <span class="hljs-number">1</span>;  
            <span class="hljs-keyword">for</span> (<span class="hljs-keyword">int</span> i = <span class="hljs-number">0</span>; i < value.length; i++) {         
                <span class="hljs-comment">//检查数组对角线位置的数值是否是0,如果是零则对该数组进行调换,查找到一行不为0的进行调换  </span>
                <span class="hljs-keyword">if</span> (value[i][i] == <span class="hljs-number">0</span>) {  
                    value = changeDeterminantNoZero(value,i,i);  
                    result*=-<span class="hljs-number">1</span>;  
                }  

                <span class="hljs-keyword">for</span> (<span class="hljs-keyword">int</span> j = <span class="hljs-number">0</span>; j <i; j++) {  
                    <span class="hljs-comment">//让开始处理的行的首位为0处理为三角形式  </span>
                    <span class="hljs-comment">//如果要处理的列为0则和自己调换一下位置,这样就省去了计算  </span>
                    <span class="hljs-keyword">if</span> (value[i][j] == <span class="hljs-number">0</span>) {  
                        <span class="hljs-keyword">continue</span>;  
                    }  
                    <span class="hljs-comment">//如果要是要处理的行是0则和上面的一行进行调换  </span>
                    <span class="hljs-keyword">if</span> (value[j][j]==<span class="hljs-number">0</span>) {  
                        <span class="hljs-keyword">double</span>[] temp = value[i];  
                        value[i] = value[i-<span class="hljs-number">1</span>];  
                        value[i-<span class="hljs-number">1</span>] = temp;  
                        result*=-<span class="hljs-number">1</span>;  
                        <span class="hljs-keyword">continue</span>;  
                    }  
                    <span class="hljs-keyword">double</span>  ratio = -(value[i][j]/value[j][j]);  
                    value[i] = addValue(value[i],value[j],ratio);  
                }  
            }           
            <span class="hljs-keyword">return</span> mathValue(value,result);
        }  

        <span class="hljs-javadoc">/** 
         * 计算行列式的结果 
         *<span class="hljs-javadoctag"> @param</span> value 
         *<span class="hljs-javadoctag"> @return</span> 
         */</span>  
        <span class="hljs-keyword">private</span> <span class="hljs-keyword">static</span> <span class="hljs-keyword">double</span> <span class="hljs-title">mathValue</span>(<span class="hljs-keyword">double</span>[][] value,<span class="hljs-keyword">double</span> result){  
            <span class="hljs-keyword">for</span> (<span class="hljs-keyword">int</span> i = <span class="hljs-number">0</span>; i < value.length; i++) {  
                <span class="hljs-comment">//如果对角线上有一个值为0则全部为0,直接返回结果  </span>
                <span class="hljs-keyword">if</span> (value[i][i]==<span class="hljs-number">0</span>) {  
                    <span class="hljs-keyword">return</span> <span class="hljs-number">0</span>;  
                }  
                result *= value[i][i];  
            }  
            <span class="hljs-keyword">return</span> result;  
        }  

        <span class="hljs-javadoc">/*** 
         * 将i行之前的每一行乘以一个系数,使得从i行的第i列之前的数字置换为0 
         *<span class="hljs-javadoctag"> @param</span> currentRow 当前要处理的行 
         *<span class="hljs-javadoctag"> @param</span> frontRow i行之前的遍历的行 
         *<span class="hljs-javadoctag"> @param</span> ratio 要乘以的系数 
         *<span class="hljs-javadoctag"> @return</span> 将i行i列之前数字置换为0后的新的行 
         */</span>  
        <span class="hljs-keyword">private</span> <span class="hljs-keyword">static</span> <span class="hljs-keyword">double</span>[] <span class="hljs-title">addValue</span>(<span class="hljs-keyword">double</span>[] currentRow,<span class="hljs-keyword">double</span>[] frontRow, <span class="hljs-keyword">double</span> ratio){  
            <span class="hljs-keyword">for</span> (<span class="hljs-keyword">int</span> i = <span class="hljs-number">0</span>; i < currentRow.length; i++) {  
                currentRow[i] += frontRow[i]*ratio;  
            }  
            <span class="hljs-keyword">return</span> currentRow;  
        }  

        <span class="hljs-javadoc">/** 
         * 指定列的位置是否为0,查找第一个不为0的位置的行进行位置调换,如果没有则返回原来的值 
         *<span class="hljs-javadoctag"> @param</span> determinant 需要处理的行列式 
         *<span class="hljs-javadoctag"> @param</span> line 要调换的行 
         *<span class="hljs-javadoctag"> @param</span> row 要判断的列 
         */</span>  
        <span class="hljs-keyword">private</span> <span class="hljs-keyword">static</span> <span class="hljs-keyword">double</span>[][] <span class="hljs-title">changeDeterminantNoZero</span>(<span class="hljs-keyword">double</span>[][] determinant,<span class="hljs-keyword">int</span> line,<span class="hljs-keyword">int</span> column){  
            <span class="hljs-keyword">for</span> (<span class="hljs-keyword">int</span> i = line; i < determinant.length; i++) {  
                <span class="hljs-comment">//进行行调换  </span>
                <span class="hljs-keyword">if</span> (determinant[i][column] != <span class="hljs-number">0</span>) {  
                    <span class="hljs-keyword">double</span>[] temp = determinant[line];  
                    determinant[line] = determinant[i];  
                    determinant[i] = temp;  
                    <span class="hljs-keyword">return</span> determinant;  
                }  
            }  
            <span class="hljs-keyword">return</span> determinant;  
        }         

        <span class="hljs-javadoc">/** 
         * 转换为代数余子式 
         *<span class="hljs-javadoctag"> @param</span> row 行 
         *<span class="hljs-javadoctag"> @param</span> line 列 
         *<span class="hljs-javadoctag"> @param</span> matrix 要转换的矩阵 
         *<span class="hljs-javadoctag"> @return</span> 转换的代数余子式 
         */</span>  
        <span class="hljs-keyword">private</span> <span class="hljs-keyword">static</span> <span class="hljs-keyword">double</span>[][] <span class="hljs-title">getNewMatrix</span>(<span class="hljs-keyword">int</span> row,<span class="hljs-keyword">int</span> line,<span class="hljs-keyword">double</span>[][] matrix){  
            <span class="hljs-keyword">double</span>[][] newMatrix = <span class="hljs-keyword">new</span> <span class="hljs-keyword">double</span>[matrix.length-<span class="hljs-number">1</span>][matrix[<span class="hljs-number">0</span>].length-<span class="hljs-number">1</span>];  
            <span class="hljs-keyword">int</span> rowNum = <span class="hljs-number">0</span>,lineNum = <span class="hljs-number">0</span>;  
            <span class="hljs-keyword">for</span> (<span class="hljs-keyword">int</span> i = <span class="hljs-number">0</span>; i < matrix.length; i++) {  
                <span class="hljs-keyword">if</span> (i == row){  
                    <span class="hljs-keyword">continue</span>;  
                }  
                <span class="hljs-keyword">for</span> (<span class="hljs-keyword">int</span> j = <span class="hljs-number">0</span>; j < matrix[i].length; j++) {  
                    <span class="hljs-keyword">if</span> (j == line) {  
                        <span class="hljs-keyword">continue</span>;  
                    }  
                    newMatrix[rowNum][lineNum++%(matrix[<span class="hljs-number">0</span>].length-<span class="hljs-number">1</span>)] = matrix[i][j];  
                }  
                rowNum++;  
            }  
            <span class="hljs-keyword">return</span> newMatrix;  
        }  

        <span class="hljs-keyword">public</span> <span class="hljs-keyword">static</span> <span class="hljs-keyword">void</span> <span class="hljs-title">main</span>(String[] args) {  
            <span class="hljs-comment">//double[][] test = {{0,0,0,1,2},{0,0,0,2,3},{1,1,0,0,0},{0,1,1,0,0},{0,0,1,0,0}};  </span>
            <span class="hljs-keyword">double</span>[][] test = {
                    {<span class="hljs-number">3.8067488033632655</span>, -<span class="hljs-number">2.894113667134647</span>},  
                    {-<span class="hljs-number">2.894113667134647</span>, <span class="hljs-number">3.6978894069779504</span>}
                };
            <span class="hljs-keyword">double</span> result;  
            <span class="hljs-keyword">try</span> {  
                <span class="hljs-keyword">double</span>[][] temp = <span class="hljs-keyword">new</span> <span class="hljs-keyword">double</span>[test.length][test[<span class="hljs-number">0</span>].length];  
                <span class="hljs-keyword">for</span> (<span class="hljs-keyword">int</span> i = <span class="hljs-number">0</span>; i < test.length; i++) {  
                    <span class="hljs-keyword">for</span> (<span class="hljs-keyword">int</span> j = <span class="hljs-number">0</span>; j < test[i].length; j++) {  
                        temp[i][j] = test[i][j];  
                    }  
                }  
                <span class="hljs-comment">//先计算矩阵的行列式的值是否等于0,如果不等于0则该矩阵是可逆的  </span>
                result = mathDeterminantCalculation(temp);  
                <span class="hljs-keyword">if</span> (result == <span class="hljs-number">0</span>) {  
                    System.out.println(<span class="hljs-string">"矩阵不可逆"</span>);  
                }<span class="hljs-keyword">else</span> {  
                    System.out.println(<span class="hljs-string">"矩阵可逆"</span>);  
                    <span class="hljs-comment">//求出逆矩阵  </span>
                    <span class="hljs-keyword">double</span>[][] result11 = getInverseMatrix(test,result);  
                    <span class="hljs-comment">//打印逆矩阵  </span>
                    <span class="hljs-keyword">for</span> (<span class="hljs-keyword">int</span> i = <span class="hljs-number">0</span>; i < result11.length; i++) {  
                        <span class="hljs-keyword">for</span> (<span class="hljs-keyword">int</span> j = <span class="hljs-number">0</span>; j < result11[i].length; j++) {  
                            System.out.print(result11[i][j]+<span class="hljs-string">"   "</span>);                       
                        }  
                        System.out.println();  
                    }  
                }  
            } <span class="hljs-keyword">catch</span> (Exception e) {  
                e.printStackTrace();  
                System.out.println(<span class="hljs-string">"不是正确的行列式!!"</span>);  
            }  
        }   
}</code>
本文内容由网友自发贡献,版权归原作者所有,本站不承担相应法律责任。如您发现有涉嫌抄袭侵权的内容,请联系:hwhale#tublm.com(使用前将#替换为@)

最小二乘法–高斯牛顿迭代法 的相关文章

  • 两直线垂直,斜率乘积为-1的证明

    老早以前在学习初等函数的时候 线性函数中的两直线y m0x b0 y m1x b1如果垂直 则有结论两条直线的斜率乘积为 1即m0 m1 1 以前也只是拿来用 没有证明过 最近在学图形学的时候 突然想起了这个点 因此记一篇笔记 证明一下 如
  • linux笔记之初次接触信号

    一 关于信号概念 1 信号是Linux所使用的进程间通信的最古老的方式 它是在软件层次上对中断机制的一种模拟 是一种异步通信的方式 一个完整的信号周期包括三个部分 信号的产生 信号在进程中的注册 信号在进程中的注销 执行信号处理函数 如下图
  • 拉格朗日插值

    直接上公式 简单的讲 这个玩意就是在给你若干个 f xi yi 的结果 算出f k 的结果 最朴素的实现方法 验证下这个公式的结果 include
  • 牛顿迭代法原理讲解

    牛顿迭代法原理讲解 牛顿迭代法是用于求解等式方程的一种方法 类似于求解F x 0的根 牛顿迭代法求解的是近似根 这个想法准确来说来源于泰勒展开式 我们知道 有些时候 我们需要求解的表达式可能非常复杂 通过一般的方法 我们很难求出它的解 所以
  • 回调函数&&回调机制

    所谓回调 定义是 一个方法的指针传递给事件源 当某一事件发生时用来调用这个方法 比如客户程序C调用服务程序S中的某个函数A 然后S又在某个时候反过来调用C中的某个函数B 对于C来说 这个B便叫做回调函数 例如Win32下的窗口过程函数就是一
  • 迭代法求解线性方程组(C++实现)

    本系列是数值分析相关算法的文章 这次用迭代法求解线性方程组 不同于上次用高斯消元法之类的求解 迭代法对于稀疏矩阵方程组的运算 会大大提高 而如果用高斯相关的算法求解 会浪费大量资源计算无用的东西 所以有必要研究此算法 本文章主要使用了3个算
  • c++11中四种类型转换

    1 static cast 功能 完成编译器认可的隐式类型转换 格式type1 a type2 b staic cast
  • 线性代数 - 特征向量和特征值

    今天在看到这个马汉诺拉距离的时候 又看到了这个东西 就是利用特征值来进行协方差方向上的伸缩 突然感觉到了线性代数的作用了 但是实际上 我今天看到了非常多的内容 但是都没有吸收完 很多内容都是线性代数的东西 但是这些东西我都忘了 这里先挖个坑
  • 三角函数常见基本公式

    定义式 图形 正弦 sin 余弦 cos 正切 tan或tg 余切 cot或ctg 正割 sec 余割 csc 函数关系 商数关系 倒数关系 平方关系 和差角公式 二角和差公式 三角和公式 积化和差公式 倍角公式 二倍角公式 三倍角公式 四
  • (邱维声)高等代数课程笔记:极大线性无关组,向量组的秩

    极大线性无关组 向量组的秩 quad 一般地 设 V V V 是数域 K K K 上的一个线性空间
  • 06. 计数原理

    6 计数原理 6 1 分类加法计数原理与分步乘法计数原理 分类加法计数原理定义 完成一件事 有 n n n 类办法 在第1类办法中有 m 1 m 1
  • 【python】numpy.percentile()函数

    numpy percentile 1 函数 百分位数是统计中使用的度量 表示小于这个值的观察值的百分比 函数numpy percentile 接受以下参数 np percentile a q axis None out None overw
  • 为什么函数y=f(x)的导数dy/dx可以适用分数运算呢?

    一 问题背景 在同济大学高等数学关于导数的内容中 如果函数y f x 可以由参数方程 表示 且三个函数皆可导 且x的值不为0 则 才开始看这个公式推导时 觉得没什么问题 仔细一想 dy dx是导数的表示符号 为什么这个符号可以适用分数运算公
  • 【华为OD机试真题 python】数字加减游戏【2022 Q4

    题目描述 数字加减游戏 小明在玩一个数字加减游戏 只使用加法或者减法 将一个数字s变成数字t 在每个回合中 小明可以用当前的数字加上或减去一个数字 现在有两种数字可以用来加减 分别为a b a b 其中b没有使用次数限制 请问小明最少可以用
  • gym 101512 BAPC 2014 I Interesting Integers

    Problem codeforces com gym 101512 attachments vjudge net contest 186506 problem I Meaning 给出一个 正整数 n 要找尽量小的 a 和 b a lt b
  • 《我的世界》Python编程入门(9) 使用函数建造房子

    一 函数的基本概念 1 1 函数在数学中的概念 函数指一个量随着另一个量的变化而变化 函数的数学形式 y f x f是一种定义好的关系 可以简称为函数 在函数f中 只要x值的确定 那么y的值一定是确定的 y的值随x值的变化而变化 1 2 P
  • 我的百度经验目录

    百度经验目录 进一步了解基于Mathematica的图像特征检测方法 http jingyan baidu com article a501d80c44a372ec630f5eb4 html 怎么把python代码打包成exe文件 http
  • Function overloaded in C++

    重载定义 如果两个函数名字相同并且在相同的域中被声明 但是参数表不同 则它们就是重载函数 重载函数条件 2 1 参数类型或参数个数不同 2 2 返回值不同不能视为重载 2 3 是否为常函数不能视为重载 2 4 对于普通类型参数只有const
  • Mathematica函数大全

    一 运算符及特殊符号 Line1 执行Line 不显示结果 Line1 line2 顺次执行Line1 2 并显示结果 name 关于系统变量name 的信息 name 关于系统变量name 的全部信息 command 执行Dos 命令 n
  • Gauss_Seidel method with python

    Gauss Seidel method with python from wikipedia https en wikipedia org wiki Gauss E2 80 93Seidel method import numpy as n

随机推荐

  • mac os x excel 单元格换行

    参考 http jingyan baidu com article 0f5fb09911cb366d8334ea07 html Windows 下是 alt 回车 mac os x 下变成 alt ctrl 回车
  • Vue3通透教程【十四】TS其他类型详解(一)

    文章目录 写在前面 对象类型 函数结构类型 数组类型 元组 枚举 类型别名 写在最后 写在前面 专栏介绍 凉哥作为 Vue 的忠实 粉丝输出过大量的 Vue 文章 应粉丝要求开始更新 Vue3 的相关技术文章 Vue 框架目前的地位大家应该
  • 毕业设计 STM32的智能WIFI视频人脸追踪监控系统

    0 前言 这两年开始毕业设计和毕业答辩的要求和难度不断提升 传统的毕设题目缺少创新和亮点 往往达不到毕业答辩的要求 这两年不断有学弟学妹告诉学长自己做的项目系统达不到老师的要求 为了大家能够顺利以及最少的精力通过毕设 学长分享优质毕业设计项
  • 函数内变量的提升(function hoisting)

    1 函数内外有重名的变量时 局部变量会覆盖全局变量 原因是函数域优先于全局域 2 当js执行进入函数时 函数内部声明过的所有变量会被提到最前 但同时对变量的赋值等操作不会被提升 例 var a 123 function test alert
  • 12帧跑步动画分解图_今天给大家分享一个跑步动画教程和注意事项!希望有所帮助!...

    跑步的动画的制作 一 跑步的基本原理 前面介绍了走路的动画的制作 跑步的制作方式和走路的方式是一样的 但是我们怎样来区别这两个动作的不同呢 虽然跑步在日常生活中经常看见 但是我们可能从来没有仔细的分析每一个动作 现在我们再来简单的说一下走路
  • upload labs第二关

    从上往下 首先定义两个变量 其中一个为空 在点击提交按钮后 前提文件路径可以找到 开始看文件类型是否为jpeg png gif格式 is upload false msg null if isset POST submit if file
  • Docker搭建zookeeper

    问题背景 前言 本文参考自 docker compose快速搭建Zookeeper集群 熬到凌晨三点多验证部署成功 网上有很多文章已经无法正确部署了 因为有些东西版本升级了 版本跟不上就会报错 还有一种更加详细更加全面的部署方式 Docke
  • 新人如何快速高效的学习Java?

    如果是新人 不想通过培训班 想学java 那么我可以很认真的告诉你 如果你是因为兴趣学学 那么你怎么学都可以 建议你找一些零基础入门的视频来学习 先看一遍 认识一下Java是个什么东西 如果是想转行学习 靠这个来工作 那么你就要好好的制定一
  • 一台计算机要两个内网,局域网如何在一台电脑上设置两个IP地址

    由于工作原因 有时需要连接两个局域网 除了频繁地更换不同局域网的网线 还要不停地设置不同局域网的IP地址 真是很麻烦 下面是学习啦小编收集整理的局域网如何在一台电脑上设置两个IP地址 希望对大家有帮助 局域网在一台电脑上设置两个IP地址的方
  • STM32F4单片机ADC采样及ARM-DSP库的FFT

    模拟信号经过ADC采样后变成数字信号 数字信号可以进行FFT运算 在频域中更容易分析信号的特征 本文将介绍如何用STM32F4的进行ADC采样 并利用ARMDSP库里的FFT算法对ADC采样值进行快速傅里叶变换 我使用的是STM32F407
  • CUDA编程中内存管理机制

    GPU设备端存储器的主要分类和特点 大小 全局 Global 和纹理 Texture 内存 大小受RAM大小的限制 本地 local 内存 每个线程限制在16KB
  • windows平台中使用vscode远程连接linux进行c++开发配置教程(内容详细适合小白)-2021-3-30

    文章目录 一 简要介绍 二 软件安装步骤 1 linux系统安装 2 vscode安装 3 ssh安装 4 配置Remote SSH 5 安装远程插件 6 简单小测试 三 配置vscode开发环境 1 默认设置 用户设置 远程设置和工作区设
  • Qt 开发环境搭建

    Qt开发环境搭建 Qt下载 Qt安装 Windows平台 离线安装 在线安装 Qt安装目录 VS2019搭建Qt开发环境 安装扩展插件 Qt VS Tools Qt Versions配置 问题 VS2019双击编辑UI时闪退 qt显示中文乱
  • 区块链物品溯源系统

    文章目录 前言 一 区块链有哪些特点 二 区块链能给品牌或者行业带来什么 1 信任度 2 小程序展示 总结 前言 区块链是一个典型的分布式协同系统 多方共同维护一个不断增长的分布式数据记录 这些数据通过密码学技术保护内容和时序 使得任何一方
  • Qt multiple definition of (function)

    前景 做项目代码优化 将原来的代码按简单工厂模式进行重新组合编写 对整个模块的文件夹进行分类 归纳 中途 出现这一问题 问题详述 某一类中的全部函数都有error multiple definition of function name 解
  • Linux 下刷 TWRP

    安装 adb 和 fastboot apt install android tools adb android tools fastboot 下载需要的 TWRP https dl twrp me flo 开机状态下进入 bootloade
  • async_await用法

    async作为修饰关键字修饰在函数前 表示该函数是一个异步函数 await的使用必须有async关键字 await等待的必须是一个promise对象 async返回的是一个promise对象 asyn function A return 星
  • pthread_cond_destroy()函数的使用

    NAME pthread cond destroy pthread cond init destroy and initialize condition variables SYNOPSIS THR include
  • 像数组一样使用NodeList:一个对象组合的有效用法

    场景 我是用querySelectorAll 查询了一些标记 并收到了一个NodeList响应 问题 节点列表类似于数组 比如 他们都有一个长度属性 它们都在括号中的索引访问它们的属性或者子元素 NodeList 0 尝试使用 map fi
  • 最小二乘法–高斯牛顿迭代法

    最小二乘法 高斯牛顿迭代法 本文将详解最小二乘法的非线性拟合 高斯牛顿迭代法 1 原理 高斯 牛顿迭代法的基本思想是使用泰勒级数展开式去近似地代替非线性回归模型 然后通过多次迭代 多次修正回归系数 使回归系数不断逼近非线性回归模型的最佳回归