数值计算实验报告

问题一：

1. 问题描述
   

2. 算法分析

   • 高斯消去法
     高斯消去法的过程与我们平常进行人工解方程组的过程相似，即通过行倍加，倍乘变换，行交换等将矩阵化简为一个三角矩阵再进行求解。化简过程通过每次迭代计算出乘数因子进行化简。以第一次迭代为例：
     第一次迭代的乘数因子为

     mi1=a(1)i1/a(1)11$$m_{i1}=a_{i1}^{(1)}/a_{11}^{(1)}$$
     
     再将 -m_{ij} 乘第一个方程，加到第i个方程。
     这就是我们常说的行倍加变换。
     需要注意的是，每次迭代的每行首位都不能为0，因此常常需要交换行

   • 列主元消去法
     列主元消去法和高斯消去法特别相似，只是增加了一些交换行的步骤。每次选取列主元中最大的一个，判断它是否为零，若非零，与第一行（相对于子矩阵）进行交换。

3. 数值实验

   1. 高斯消去法
      编写函数guasselim进行高斯消元，并计算其运行时间。迭代过程通过循环完成

          for k =1:n-1;
              %compute the kth column of M
              m(k+1:n) = A(k+1:n,k)/A(k,k);
              %compute
              %An=Mn*An-1
              %bn=Mn*bn-1;
              for i=k+1:n
                  A(i,k+1:n) = A(i,k+1:n)-m(i)*A(k,k+1:n);
              end;
              b(k+1:n) = b(k+1:n)-b(k)*m(k+1:n);
          end

   2. 列主元消去法
      列主元消去法的特点在于对最大列主元的选择以及行交换，相当于对高斯消元法的一种改进。由于matlab提供了十分方便的寻找最大元素的函数，以及十分简便的交换语法，因此执行上述操作非常方便

           %find the max element of the kth column
           [max_a,index] = max(abs(A(1:n,k))); %取第k列绝对值最大数的值和索引
           if max_a == 0
               mdet = 0;
               return;
           end;
           A([k,index],:) = A([index,k],:);
           b([k,index],:) = b([index,k],:);  %交换k行和index行
           mdet = -mdet;

4. 结果测试分析
   编写了测试函数，对满足正态分布的矩阵和向量A，b分别采用高斯消元法和列主元消去法进行求解，矩阵大小n = [10,50,100,200].
   结果如下：
   
   多次实验发现（这里不一一列出），列主元消去法的时间通常要比高斯消元法的时间要多，原因可能在于，列主元需要时时进行行交换，导致效率有点下降。不过，与高斯消元法相比较，列主元消去法的结果要比高斯消元法的结果更好。因为列主元消去法每次都会将最大的第一位元素所在的行移到第一行，于是计算乘数因子的时候，不会出现分母过小导致的误差变大的问题

---

问题二：

1. 问题描述：
   

2. 算法分析
   说到迭代求解线性方程组，首先要讲的是一个基本方法。此后的雅各比迭代，高斯塞德尔迭代，超松弛迭代，都是基于这个方法进行操作的，
   假设有线性方程组
   

   Ax=b$$Ax=b$$
   
   将A分裂为
   
   A=M−N$$A=M-N$$
   
   则原线性方程组可以转化为求解
   
   Mx=Nx+b$$Mx=Nx+b$$
   
   即
   
   x=M−1Nx+M−1b$$x=M^{-1}Nx+M^{-1}b$$
   
   根据上式，我们就能够进行一阶定常迭代

   接下来，根据不同的算法来做具体分析

   1. 雅各比迭代法
      雅各比迭代法将A 化为 D - L - U， D是对角矩阵， L 和 U分别是下三角矩阵和上三角矩阵
      取M = D， N = L + U， 满足 A = M - N
      利用

      x=M−1Nx+M−1b$$x=M^{-1}Nx+M^{-1}b$$
      进行迭代计算

   2. 高斯塞德尔迭代法
      令M = D - L， N = U， 同上进行迭代运算

   3. 超松驰迭代法
      

      M=1w(D−wL)$$M=\frac{1}{w}(D-wL)$$
      , w称为松弛因子。代入迭代式中就能进行迭代。需要注意的是，松弛银子的范围设为（0，2）。当松弛因子值为1时，SOR方法即为高斯塞戴尔迭代法

   4. 共轭梯度法
      CG方法是一种十分常用的变分方法。所谓变分法，就是为了最终寻求极值函数，对应于求一个二次函数的极值。
      CG方法用于求解大型稀疏对称矩阵非常有效。

3. 数值实验

   1. 雅各比迭代法

      for j=1:N
          x(j) = (b(j)-A(j,[1:j-1,j+1:N]) * P([1:j-1, j+1:N]))/A(j,j);
      end

   2. 高斯塞德尔迭代法

          for i=1:n
      
              sigma=0;
      
              for j=1:i-1
                      sigma=sigma+A(i,j)*x(j);
              end
      
              for j=i+1:n
                      sigma=sigma+A(i,j)*x_old(j);
              end
      
              x(i)=(1/A(i,i))*(b(i)-sigma);
          end

   3. 超松驰迭代法

      D = diag(diag(A));
      L = -tril(A,-1);
      U = -triu(A,1);
      B = (D - L * w)\((1 - w) * D + w * U);
      f = w * inv((D - L * w)) * b;
      x = B * x0 + f;
      n = 1;

   4. 共轭梯度法

          for k = 1:numel(b);
             r = r - t*z;
             err = norm(r);
             if( norm(r) < tol )
                  return;
             end
             B = (r'*z)/s;
             y = -r + B*y;
             z = A*y;
             s = y'*z;
             t = (r'*y)/s;
             x = x + t*y;
          end

4. 结果测试分析。
   
   
   由结果来看，随着n的增加，Jacobi 迭代法的运行效率持续下降，甚至从收敛变为发散。高斯塞戴尔次之，而共轭梯度法表现最好。对比上下两图，区别在于我将y轴范围扩大。可以看出，雅各比迭代法在迭代早期就开始发散。 而观察共轭梯度法，当n较小时，共轭梯度法只需要几步迭代就能接近最优解。而随着n增加，仍然保持着较好的效果。原因在于我们将A设为对称正定矩阵，而共轭梯度法十分适合这类矩阵

   观察以下时间运行比较
   
   可见，共轭梯度法仍然是最优秀的算法

ps.如何生成一个特征值在0到1之间均匀分布的对称正定矩阵？

function [A] = generateSPDmatrix(n)
%A = rand(n,n); 
%A = 0.5*(A+A'); 
%A = A + n*eye(n); 
d = rand(1,n);
D = diag(d); %d作为对角元素
P = rand(n);  
[Q,~] = qr(P); %将P进行QR分解获得正交矩阵
A = Q*D*Q';  %对对角矩阵D进行相似变换获得A

end

---

问题三

1. 问题描述
   

   给定的网页显示如下信息
   

2. 算法分析

   1. 基本概念
      pagerank算法也是通过迭代来对多个节点进行排名。假设P为结果向量，向量元素为P的价值。那么有
      

      Pn+1=APn$$P_{n+1}=A{P}_n$$
      
      
      A=dS+1−aNeeT$$A=dS+\frac{1-a}{N}e{e}^T$$
      • 矩阵S记录了不同节点之间的关系，从一个例子中来说明：
        
        引用自pagerank从原理到实现
      • a表示权重，一般为0.85
      • (1-a) / N是为没有出度的节点设置的，假设他们对其他所有节点都有贡献（模拟了用户进入某个不链接到其他网页的网页上时，会通过输入网址随机访问其他网页）

   2. 实现过程
      显然，基本概念中的描述是对矩阵进行操作，可是本次实验要求给出的节点数为75879个，边数为508837。当你尝试着构造S矩阵时你会发现，matlab报错，原因提示为75879 * 75879 矩阵需要占据50G的内存，超出最大限制。这种方法行不通。后来与他人交流，发现能够使用稀疏矩阵进行操作，不过下面要讲的这种方法将更加简便。
      从上述的S矩阵例子可以看出，矩阵的每一列凡是非零的元素，值都相等。原因在于，sij表示的是j对于i的出链，这里假设为出链之和是1，出链的概率相同。因此，假如j有n个出链，那么每个出链的贡献即为1/n
      因此，我们可以得出我们的实现过程。

      1. 读入数据文件，得到一个75879 * 2 的矩阵
      2. 遍历数据矩阵的第一列，计算每种元素的个数并存入一个75879长度的数组。所获得数组即为每种元素的出度数量。根据出度数量就可以获得每种元素每条出链的贡献值
      3. 遍历数据矩阵的第二列，对于不同的数字j，我们可以找到对应的第一列数字i，这个数字即为i对j的入度。根据i去上一步的数组中寻找对应的贡献值，与旧结果数组的对应位置相乘后加到新的结果数组的对应元素中。
      4. 最后一步，即为结果数组中的每一个值添加一个 （1-a）/ N * x 。x为对应旧结果数组的每个值。

      可以发现，我们的实现方法完全模拟矩阵与向量相乘，但是这是在不构造矩阵的情况下完成。因为发现（1-a）/ N 对于任何元素（不管有没有出度），因此留到了最后直接与所有旧数组元素相乘之后加到新数组的每个元素。这个实现过程很简单，调用sum函数计算旧数组中所有元素之和，再乘以（1-a） / N

   3. 实现过程中的问题。

      1. 寻找每个元素入度过程中出现的问题
         这是未修改前的代码：

         %    for i = 1:MAXINDEX  %对结果数组的每一个进行操作
         %        index = i - 1;
         %        tempResult = 0;
                 for j = 1:EDGE
                     index = soc_Epinions1(j,2);
                     target = soc_Epinions1(j,1);  %target是从0开始
                     tempA = a * (1/ODnum(target + 1));
                     tempJ = tempA * oldResult(target + 1);
                     Result(index+1) = Result(index+1) + tempJ;
                 end
         %        tempResult = tempResult + ((1-a)/N) * sum(oldResult);
         %        Result(i) = tempResult;
                 for j = 1:MAXINDEX
                     if Result(j) ~= 0
                         Result(j) = Result(j) + ((1-a)/N) * sum(oldResult);
                     end
                 end
         
                 change = change + abs(sum(Result - oldResult));
         %    end

         从代码可以看出，为了寻求不同元素的入度，我在外层（变量i所在层）进行循环，每次选择特定的元素i，然后遍历数据矩阵去寻找i的入度。数据矩阵有50万条数据，节点有7万个，因此一次迭代下来，我需要做350万次计算
         毫无疑问，结果是跑不出来的
         因此我将代码做出以下修改：

                 for j = 1:EDGE
                     index = soc_Epinions1(j,2);
                     target = soc_Epinions1(j,1);  %target是从0开始
                     tempA = a * (1/ODnum(target + 1));
                     tempJ = tempA * oldResult(target + 1);
                     Result(index+1) = Result(index+1) + tempJ;
                 end
                 for j = 1:MAXINDEX
                     if Result(j) ~= 0
                         Result(j) = Result(j) + ((1-a)/N) * sum(oldResult);
                     end
                 end

         思想就在于，不再选择特定的元素，而是在扫描数据矩阵的第二列元素的时候，扫到什么元素就算什么元素，将计算结果加到Result对应元素里面。因此，一次扫描就能完成所有元素的值的更新。最后，再将 (1-a) / N 加到Result的每个元素里。这里有需要一次循环。因此总共迭代大概为57万次。

      2. 数据矩阵的下标是从0开始的，而matlab数组的下标从1开始。
         这个问题，只需要用进行1的增加或减少即可。在代码中，我用index表示数据下标（从0开始），用i表示数组下标（从1开始）

      3. 数据节点数为75879个，而查询数据文件发现，最大下标为75887。这是，我开始担心，节点标号并不是连续的。毫无疑问，我必须申请一个长度大于75888的数组来存储，但因此也必然会出现，数组中的某些元素，在实际数据文件中并不存在，因此迭代过程中必须让这些元素保持为0（就如同不存在这个网页，就不能对他进行操作）。因此，在初始化和迭代过程中，我小心的避免了这些节点

         for i = 1:EDGE  %初始化
             index1 = soc_Epinions1(i,1); %index 从零开始
             index2 = soc_Epinions1(i,2);
             Result(index1 + 1) = 1/N;
             Result(index2 + 1) = 1/N;
             ODnum(index1 + 1) = ODnum(index1 + 1) + 1;
         end

         
                 for j = 1:MAXINDEX  %迭代
                     if Result(j) ~= 0
                         Result(j) = Result(j) + ((1-a)/N) * sum(oldResult);
                     end
                 end

   4. 结果分析
      
      由下至上，排名逐渐降低，第一列为元素标号，第二列为对应价值。从图中看出，第一名为18号元素
      精度设置为 1e-4
      迭代次数为1次