跳转至

基础数值算法

数值问题与误差

数值问题的特点

数值算法处理连续量,通常求的是满足误差要求的近似解。除了“有没有解、解是否唯一”,还要问:输入稍有变化时解会不会剧烈变化,有限精度计算是否稳定,以及达到目标精度需要多少计算量。

Hadamard 意义下的适定问题满足三个条件:解存在、解唯一、解连续依赖输入。缺少任何一项就是不适定问题。适定仍不等于良态;“连续依赖”只排除了跳变,误差放大倍数仍可能很大。

问题的条件性与算法的稳定性要分开:

  • 病态问题:输入的微小扰动就可能引起很大的解变化,这是问题本身的性质;
  • 不稳定算法:计算过程中又额外放大了舍入、截断等误差;
  • 稳定算法只能避免雪上加霜,不能把病态问题变成良态问题。

误差来源与表示

设真值为 \(x\)、近似值为 \(\tilde x\),绝对误差和相对误差为

\[ e=\tilde x-x, \qquad e_r=\frac{\tilde x-x}{x}. \]

讨论误差大小时使用 \(|e|,|e_r|\)。当真值接近 \(0\) 时,相对误差失去意义,应改看绝对误差或按问题尺度归一化。

数值误差有几种不同来源:

  • 模型误差:数学模型省略了现实因素;
  • 测量/数据误差:输入本身不精确;
  • 截断误差:用有限项级数、有限差分或有限次迭代代替无限过程;
  • 舍入误差:实数只能用有限位浮点数表示,每次运算还要再次舍入。

对精确映射 \(f\)、近似输入 \(\tilde x\) 和实际计算过程 \(\tilde f\),总误差可拆为

\[ f(x)-\tilde f(\tilde x) =\underbrace{f(x)-f(\tilde x)}_{\text{输入误差的传播}} +\underbrace{f(\tilde x)-\tilde f(\tilde x)}_{\text{算法与舍入误差}}. \]

有限差分:步长不是越小越好

Taylor 展开给出前向差分与中心差分:

\[ f'(x)=\frac{f(x+h)-f(x)}{h}+O(h), \]
\[ f'(x)=\frac{f(x+h)-f(x-h)}{2h}+O(h^2). \]

若邻域内 \(|f''|\le M\),函数值的绝对误差不超过 \(\varepsilon\),则前向差分的总误差可以估成

\[ E(h)\le \frac{Mh}{2}+\frac{2\varepsilon}{h}. \]

第一项随 \(h\) 变小而下降,第二项却因相近数相减并除以小量而上升。误差界在

\[ h\approx2\sqrt{\frac{\varepsilon}{M}} \]

附近最小。这是数值计算里很常见的权衡:离散得过粗有截断误差,取得过细又会让舍入误差占主导。

前向误差、后向误差与条件数

计算 \(y=f(x)\) 得到 \(\hat y\)。前向误差直接比较 \(\hat y\)\(y\);后向误差则问:能否找到一个邻近输入 \(\hat x=x+\Delta x\),使得 \(f(\hat x)=\hat y\)

标量函数在 \(x\) 附近的绝对、相对条件数分别为

\[ \kappa_{\rm abs}(x)=|f'(x)|, \qquad \kappa_{\rm rel}(x)=\left|\frac{x f'(x)}{f(x)}\right|. \]

一阶近似下

\[ \text{相对前向误差} \lesssim \kappa_{\rm rel}(x)\times\text{相对后向误差}. \]

后向稳定算法给出的结果,可以解释成某个非常接近的输入的精确答案。只有“问题良态 + 算法稳定”同时成立,才有理由期待较小的前向误差。

浮点数系统

抽象浮点系统 \(F(\beta,p,L,U)\) 中,正规数写成

\[ x=\pm\left(d_0+\frac{d_1}{\beta}+\cdots+ \frac{d_{p-1}}{\beta^{p-1}}\right)\beta^E, \quad \beta\ge2,\quad d_0\ne0,\quad d_i\in\{0,1,\ldots,\beta-1\},\quad L\le E\le U. \]

基数为 \(\beta\),尾数精度为 \(p\) 位,指数范围决定上溢和下溢边界。次正规数允许最小指数下的首位为零,用逐渐降低精度换取“逐渐下溢”。

IEEE 754 单精度由 \(1\) 位符号、\(8\) 位指数、\(23\) 位显式小数组成,实际有效精度为 \(24\) 个二进制位;双精度为 \(1+11+52\),有效精度为 \(53\) 位。特殊编码还表示 \(\pm0\)\(\pm\infty\) 和 NaN。

常用舍入规则有两种:向零截断直接丢弃第 \(p\) 位以后的数字;最近舍入取距离最近的浮点数,恰好居中时令末位为偶数,这是 IEEE 754 的默认方式。相应机器精度为

\[ \varepsilon_{\rm mach}= \begin{cases} \beta^{1-p},&\text{向零截断},\\ \frac12\beta^{1-p},&\text{最近舍入}. \end{cases} \]

二进制单、双精度分别约为 \(2^{-24}\approx10^{-7}\)\(2^{-53}\approx10^{-16}\)。正常范围内可用模型

\[ \operatorname{fl}(x)=x(1+\delta), \qquad |\delta|\le\varepsilon_{\rm mach}. \]

\(x=\pm m\beta^E\),近似值 \(x^*\) 满足

\[ |x-x^*|\le\frac12\beta^{E-n+1}, \]

就称它从最左边起有 \(n\) 位有效数字。若首位有效数为 \(d_0\),相对误差有

\[ \frac{|x-x^*|}{|x|} \le\frac{1}{2d_0}\beta^{-(n-1)}. \]

浮点集合对四则运算不封闭,结合律、分配律也可能失效。两个接近的大数相减会损失大量有效数字,称为灾难性抵消。二次方程可改写为

\[ q=-\frac12\left(b+\operatorname{sign}(b)\sqrt{b^2-4ac}\right), \qquad x_1=\frac qa,\quad x_2=\frac cq, \]

这里约定 \(\operatorname{sign}(0)=1\)。该写法避免在 \(-b\pm\sqrt\Delta\) 中直接相减。大量不同量级的正数求和时,通常从小到大累加,或采用补偿求和,能减少“小数被大数吃掉”的现象。

线性方程组

解、范数与条件数

一般形式是

\[ A x=b,\qquad A\in\mathbb R^{m\times n},\quad b\in\mathbb R^m. \]

方程相容当且仅当 \(b\) 属于 \(A\) 的列空间。对方阵,\(A\) 非奇异时有唯一解;\(A\) 奇异且相容时有无穷多解,奇异且不相容时无解。欠定或超定只描述方程数和未知量数的关系,不能单独决定解的个数。

常用向量范数为

\[ \|x\|_1=\sum_i|x_i|, \qquad \|x\|_2=\sqrt{\sum_i|x_i|^2}, \qquad \|x\|_\infty=\max_i|x_i|. \]

矩阵的诱导范数定义为

\[ \|A\|_p=\max_{x\ne0}\frac{\|Ax\|_p}{\|x\|_p}. \]

其中 \(\|A\|_1\) 是最大绝对列和,\(\|A\|_\infty\) 是最大绝对行和,

\[ \|A\|_2=\sigma_{\max}(A)=\sqrt{\lambda_{\max}(A^TA)}. \]

Frobenius 范数则为 \(\|A\|_F=\sqrt{\sum_{ij}|a_{ij}|^2}\),它不是同一向量范数诱导出来的算子范数。

非奇异方阵的条件数为

\[ \kappa(A)=\|A\|\,\|A^{-1}\|, \qquad \kappa_2(A)=\frac{\sigma_{\max}(A)}{\sigma_{\min}(A)}. \]

\(\kappa(A)\ge1\);奇异矩阵的条件数视为无穷。只扰动右端时,一阶误差满足

\[ \frac{\|\Delta x\|}{\|x\|} \le \kappa(A)\frac{\|\Delta b\|}{\|b\|}. \]

计算解 \(\hat x\) 的残差是 \(r=b-A\hat x\)。由 \(\hat x-x=-A^{-1}r\) 可得

\[ \frac{\|\hat x-x\|}{\|\hat x\|} \le \kappa(A)\frac{\|r\|}{\|A\|\,\|\hat x\|}. \]

所以小残差只说明后向误差小;矩阵病态时,解仍可能相差很大。若输入相对精度约为机器精度,通常会损失约 \(\log_{10}\kappa(A)\) 个十进制有效位。

高斯消去与选主元

\(i\) 步用主元 \(a_{ii}\) 消去其下方元素。对 \(k=i+1,\ldots,n\),令

\[ m_{ki}=\frac{a_{ki}}{a_{ii}}, \]

再更新

\[ a_{kj}\leftarrow a_{kj}-m_{ki}a_{ij}, \qquad b_k\leftarrow b_k-m_{ki}b_i. \]

得到上三角系统后,从后向前回代:

\[ x_i=\frac{b_i-\sum_{j=i+1}^{n}u_{ij}x_j}{u_{ii}}. \]

稠密消元为 \(O(n^3)\),回代为 \(O(n^2)\)。矩阵可以原地覆盖,但主元为零会除零,主元很小则会产生很大的乘子并放大舍入误差。

部分选主元在当前列尚未处理的行中选择绝对值最大的元素并交换到主元位置,使 \(|m_{ki}|\le1\)。完全选主元在整个剩余子矩阵中选最大元素,需要同时交换行和列;列交换改变未知量次序,最后必须还原。系数跨度很大时可先作行列缩放,但数值判秩应使用与问题尺度匹配的容差,不能机械比较是否恰等于零。

消去矩阵与 LU 分解

\(k\) 次消元可写成单位下三角消去矩阵 \(M_k\)。连续消元满足

\[ M_{n-1}\cdots M_1A=U, \]

因而不交换行时 \(A=LU\);含部分选主元时应写成

\[ PA=LU. \]

解方程分成两次三角求解:

\[ Ly=Pb,\qquad Ux=y. \]

分解仍需 \(O(n^3)\),但同一个 \(A\) 对多个右端项时只分解一次,每个新右端只花 \(O(n^2)\)。不要显式求 \(A^{-1}\) 再乘 \(b\),这比三角求解更慢,也通常更不稳定。

高斯-约当法把主元上下方都消成零,不再需要回代,但常数开销更大。若 \(A\) 对称正定,则应优先使用 Cholesky 分解

\[ A=LL^T, \]
\[ l_{jj}=\sqrt{a_{jj}-\sum_{k<j}l_{jk}^2}, \qquad l_{ij}=\frac{a_{ij}-\sum_{k<j}l_{ik}l_{jk}}{l_{jj}}. \]

它不需选主元、只存半个矩阵,主要工作量约为普通高斯消去的一半。开方量非正说明矩阵不是正定,或计算已经受到严重数值误差影响。

迭代法

\(A=M-N\),线性方程可写成不动点形式

\[ x=Gx+c,\qquad G=M^{-1}N,\quad c=M^{-1}b, \]

并迭代

\[ x^{(k+1)}=Gx^{(k)}+c. \]

对任意初值收敛的充要条件是谱半径

\[ \rho(G)=\max_i|\lambda_i(G)|<1. \]

\(A=D+L+U\),其中 \(D\) 是对角部分,\(L,U\) 是严格下、上三角部分。Jacobi 法全部使用上一轮分量:

\[ x^{(k+1)}=D^{-1}\bigl[b-(L+U)x^{(k)}\bigr], \]
\[ x_i^{(k+1)}= \frac{b_i-\sum_{j\ne i}a_{ij}x_j^{(k)}}{a_{ii}}. \]

Gauss-Seidel 法立即使用本轮已经算出的新分量:

\[ (D+L)x^{(k+1)}=b-Ux^{(k)}, \]
\[ x_i^{(k+1)}= \frac{b_i-\sum_{j<i}a_{ij}x_j^{(k+1)}- \sum_{j>i}a_{ij}x_j^{(k)}}{a_{ii}}. \]

严格行对角占优可保证两者收敛;对称正定时 Gauss-Seidel 收敛。每轮稠密计算为 \(O(n^2)\),稀疏矩阵则约为 \(O(\operatorname{nnz}(A))\)。SOR、最速下降和共轭梯度是同一方向的扩展。

停止时应同时检查归一化残差、相邻迭代差和最大迭代次数。只看 \(\|x^{(k+1)}-x^{(k)}\|\),可能把数值停滞误判成已经收敛。

数据拟合与插值

最小二乘拟合

给定 \(m\) 个观测 \((t_i,y_i)\),用基函数 \(\phi_j\) 构造模型

\[ f(t,x)=\sum_{j=1}^{n}x_j\phi_j(t), \qquad A_{ij}=\phi_j(t_i),\quad b_i=y_i. \]

\(m>n\) 时通常无法让模型通过全部数据点,于是求

\[ \min_x\|b-Ax\|_2^2. \]

记残差 \(r=b-Ax\)。最优点处,残差与 \(A\) 的列空间正交:

\[ A^Tr=0, \]

从而得到正规方程

\[ A^TAx=A^Tb. \]

最小二乘问题总有解;若 \(A\) 满列秩则解唯一。满列秩时

\[ A^+=(A^TA)^{-1}A^T, \qquad x=A^+b, \]

\(P=A(A^TA)^{-1}A^T\) 是把 \(b\) 正交投影到 \(\mathcal R(A)\) 的投影矩阵。

正规方程可以先形成 \(A^TA,A^Tb\),再对对称正定矩阵作 Cholesky 分解。它速度较快,但

\[ \kappa_2(A^TA)=\kappa_2(A)^2, \]

会明显恶化条件性。也可把残差作为未知量,解增广系统

\[ \begin{bmatrix} I&A\\ A^T&0 \end{bmatrix} \begin{bmatrix}r\\x\end{bmatrix} = \begin{bmatrix}b\\0\end{bmatrix}, \]

但矩阵更大且对称不定,需要合适的主元策略。

QR 分解与 Householder 变换

\[ A=Q\begin{bmatrix}R\\0\end{bmatrix}, \qquad Q^TQ=I, \qquad Q^Tb=\begin{bmatrix}c_1\\c_2\end{bmatrix}, \]

\[ \|b-Ax\|_2^2=\|c_1-Rx\|_2^2+\|c_2\|_2^2. \]

因此只需回代解 \(Rx=c_1\),最小残差范数为 \(\|c_2\|_2\)。QR 不形成 \(A^TA\),不会把条件数平方,通常比正规方程稳健。

Householder 反射写成

\[ H=I-2\frac{vv^T}{v^Tv}, \qquad H^T=H=H^{-1}. \]

对待消元向量 \(a\),取

\[ \alpha=-\operatorname{sign}(a_1)\|a\|_2, \qquad v=a-\alpha e_1, \]

其中同样约定 \(\operatorname{sign}(0)=1\)。这样可以把 \(a\) 反射到坐标轴,并避免 \(a_1\)\(\|a\|\) 相减造成抵消。实现时不必显式形成 \(H\) 或完整的 \(Q\),只需保存反射向量并同步作用于 \(A\)\(b\)

Householder QR 的主要工作量约为 \(mn^2-n^3/3\)。缺秩或近缺秩时,可采用带列主元的 QR 判断数值秩;若还要求最小范数解,SVD 更可靠但更昂贵。

课件中的五点二次拟合为

\[ (t_i,y_i)=(-1,1),(-0.5,0.5),(0,0),(0.5,0.5),(1,2), \]

\(1,t,t^2\) 为基底,得到

\[ f(t)\approx0.0857+0.400t+1.43t^2. \]

插值与拟合的区别

插值要求 \(p(x_i)=y_i\),拟合只要求整体残差尽量小。有噪声的数据若强制逐点通过,往往会把噪声也当成规律;拟合通常更合适。插值还要关心点间的光滑性、单调性、凸性和数值稳定性,不能只检查节点上的误差为零。

用单项式基底

\[ p_{n-1}(x)=c_0+c_1x+\cdots+c_{n-1}x^{n-1} \]

得到 Vandermonde 矩阵。节点互异时

\[ \det V=\prod_{1\le i<j\le n}(x_j-x_i)\ne0, \]

故存在唯一的 \(n-1\) 次插值多项式。不过 Vandermonde 矩阵常随节点数和尺度迅速变得病态。改变基底不改变精确算术下的多项式,却会明显影响系数求解的稳定性。

多项式求值应使用 Horner 形式

\[ p(x)=c_0+x\bigl(c_1+x(c_2+\cdots+xc_{n-1})\bigr), \]

只需 \(O(n)\) 次乘加。

Lagrange 与 Newton 插值

Lagrange 基函数为

\[ \ell_j(x)=\prod_{k\ne j}\frac{x-x_k}{x_j-x_k}, \qquad \ell_j(x_i)=\delta_{ij}, \]

于是

\[ p_{n-1}(x)=\sum_{j=1}^{n}y_j\ell_j(x). \]

它形式直接,但朴素求值较慢,增加新节点时也不便复用已有计算。

Newton 形式采用逐步增加的基底:

\[ p_{n-1}(x)=c_0+c_1(x-x_0)+c_2(x-x_0)(x-x_1)+\cdots. \]

系数由差商表在 \(O(n^2)\) 时间内得到,新节点可增量加入;嵌套求值仍为 \(O(n)\)。单项式、Lagrange 和 Newton 表示在精确算术下是同一个唯一多项式,区别在构造、更新、求值和数值性质。

高阶全局多项式容易振荡、条件性也差,实际中常改用分段三次样条、B 样条、正交多项式、有理函数或径向基函数。

非线性方程

根、条件性与收敛阶

一元问题是求 \(f(x)=0\) 的实根。若 \(f\)\([a,b]\) 连续且 \(f(a)f(b)<0\),中值定理保证区间内至少有一根,但不保证唯一,也检测不到不变号的偶重根。

\[ f(x^*)=f'(x^*)=\cdots=f^{(m-1)}(x^*)=0, \qquad f^{(m)}(x^*)\ne0, \]

\(x^*\)\(m\) 重根。单根满足 \(f'(x^*)\ne0\),它对函数值扰动的绝对条件数约为

\[ \kappa_{\rm abs}\approx\frac{1}{|f'(x^*)|}. \]

切线近水平时,小残差未必意味着位置接近真根。

若误差 \(e_k=x_k-x^*\) 满足

\[ \lim_{k\to\infty}\frac{|e_{k+1}|}{|e_k|^r}=C>0, \]

则称收敛阶为 \(r\)\(r=1,C<1\) 是线性收敛,\(r>1\) 是超线性,\(r=2\) 是平方收敛。

二分法

在连续且端点异号的区间中,每次取中点并保留仍异号的一半。\(k\) 次后区间长度为

\[ \frac{b-a}{2^k}, \]

要使它不超过 tol,需

\[ k\ge\left\lceil\log_2\frac{b-a}{\mathrm{tol}}\right\rceil. \]

二分法线性收敛,每轮约增加一个二进制有效位。它慢,但只用函数值符号且始终保留根的包围区间。中点写成 \(a+(b-a)/2\) 可降低溢出风险;实现还应处理端点恰为根、NaN 和最大迭代次数。

不动点迭代

把方程改写成 \(x=g(x)\),再计算

\[ x_{k+1}=g(x_k). \]

局部误差满足 \(e_{k+1}=g'(\xi_k)e_k\)。若 \(|g'(x^*)|<1\),根附近存在收敛域;若 \(|g'(x^*)|>1\) 则局部发散。同一个方程用不同方式改写,可能得到完全不同的收敛行为。

更强的区间保证是

\[ g([a,b])\subset[a,b], \qquad \sup_{x\in[a,b]}|g'(x)|\le q<1. \]

此时从区间内任意初值出发,都收敛到唯一不动点。

Newton、割线与保护法

Newton 法用切线代替函数:

\[ x_{k+1}=x_k-\frac{f(x_k)}{f'(x_k)}. \]

对单根,若初值足够近且函数二阶光滑,通常平方收敛。初值远、导数很小或根为重根时,可能震荡、越界、转向其他根或退化为线性收敛。

不便求导时可用割线法:

\[ x_{k+1}=x_k-f(x_k) \frac{x_k-x_{k-1}}{f(x_k)-f(x_{k-1})}. \]

它需要两个初值,每轮只新增一次函数求值,收敛阶约为黄金比 \(1.618\)。分母接近零时必须保护。反二次插值用三个点拟合反函数 \(x=p(y)\),以 \(p(0)\) 作为新点,局部收敛阶约为 \(1.839\),但同样缺少全局保证。

实用求根器通常维护一个端点异号区间:Newton、割线或反插值给出的候选点若越界或进展不足,就退回二分。Brent 类方法正是这种“快速步 + 安全步”的组合。

停止条件最好同时看区间宽度、步长、残差、导数/分母是否过小和迭代上限。病态根附近尤其不能只看 \(|f(x)|\)

最优化初步

优化问题、存在性与凸性

一般问题写成

\[ \min_{x\in S} f(x), \]

带连续约束时常写作

\[ \min f(x) \quad\text{s.t.}\quad g(x)=0,\quad h(x)\le0. \]

局部最小只在某个邻域内不比周围点差,全局最小则对所有可行点成立。连续函数在非空紧集上一定取得最小值;若可行集非空、闭但无界,而 \(f\) 连续且 \(\|x\|\to\infty\) 推出 \(f(x)\to+\infty\),也可保证最小值存在。

\(S\) 为凸集,且

\[ f(\alpha x+(1-\alpha)y) \le \alpha f(x)+(1-\alpha)f(y), \quad 0\le\alpha\le1, \]

\(f\) 是凸函数。凸优化中任何局部最小都是全局最小;严格凸时全局最小至多一个。

一阶与二阶条件

无约束内点局部极小的必要条件为

\[ \nabla f(x^*)=0. \]

若二阶可导,局部极小还要求 Hessian \(H_f(x^*)\) 半正定;若梯度为零且 Hessian 正定,则是严格局部最小。Hessian 不定时是鞍点,半正定或奇异时通常还需看更高阶项。

一维最小点附近

\[ f(x^*+h)-f(x^*)\approx\frac12f''(x^*)h^2. \]

因此函数值误差只有 \(\varepsilon\) 时,位置误差可能仍有

\[ |h|\approx\sqrt{\frac{2\varepsilon}{|f''(x^*)|}}. \]

曲率小的极小点在位置上更敏感。

等式约束问题使用 Lagrange 函数

\[ L(x,\lambda)=f(x)+\lambda^Tg(x). \]

在常用的约束资格条件下,例如 \(J_g(x^*)\) 满行秩,一阶必要条件为

\[ \nabla f(x^*)+J_g(x^*)^T\lambda^*=0, \qquad g(x^*)=0. \]

二阶条件要限制在可行切空间 \(J_gd=0\) 上判断,而不是笼统要求整个 Lagrange Hessian 正定。不等式约束则进一步引出 KKT 条件。

一维搜索

\(f\)\([a,b]\) 单峰,可用黄金分割搜索。令

\[ \tau=\frac{\sqrt5-1}{2}\approx0.618, \]
\[ x_1=a+(1-\tau)(b-a), \qquad x_2=a+\tau(b-a). \]

比较 \(f(x_1),f(x_2)\) 后舍弃不可能包含极小点的一端。黄金比例保证下一轮复用一个旧点,所以每轮只新增一次函数求值;区间长度每轮乘 \(\tau\)

逐次抛物插值用三个点拟合二次函数,以拟合抛物线的极小点作为候选,局部收敛阶约为 \(1.324\),但不保证落在可信区间。实用实现常把它和黄金分割组合。若可求二阶导,则可对 \(f'(x)=0\) 使用 Newton:

\[ x_{k+1}=x_k-\frac{f'(x_k)}{f''(x_k)}. \]

多维方法

直接搜索只用函数值,不依赖导数,适合低维或不光滑问题;维数升高后,需要探索的方向和函数求值数很快增加。

最速下降沿负梯度方向:

\[ x_{k+1}=x_k-\alpha_k\nabla f(x_k), \]

步长 \(\alpha_k\) 由线搜索决定。它稳健但通常只有线性收敛,在狭长谷底会来回折线前进。更一般的下降方向 \(d_k\) 应满足 \(\nabla f(x_k)^Td_k<0\)

Newton 法每轮解线性系统

\[ H_f(x_k)s_k=-\nabla f(x_k), \qquad x_{k+1}=x_k+s_k. \]

不要显式形成 \(H_f^{-1}\)。它在非退化极小点附近平方收敛,但远离解时 Hessian 可能不正定,Newton 方向也未必下降,需要线搜索、阻尼或 Hessian 修正。

非线性最小二乘

\[ f(x)=\frac12\|g(x)\|_2^2 \]

可用 Gauss-Newton 步

\[ (J^TJ)s=-J^Tg, \]

它忽略了 Hessian 中含残差二阶导的项,适合 Jacobian 满列秩且残差不太大的情形。拟 Newton 法则用逐步更新的 \(B_k\) 近似 Hessian;BFGS 是常用代表,以每轮 \(O(n^2)\) 的更新换取较好的全局表现和局部超线性收敛。大规模问题还常用共轭梯度与预条件技术。

评论