We introduce innovative algorithms for computing exact or approximate (minimum-norm) solutions to $Ax=b$ or the {\it normal equation} $A^TAx=A^Tb$, where $A$ is an $m \times n$ real matrix of arbitrary rank. We present more efficient algorithms when $A$ is symmetric PSD. First, we introduce the {\it Triangle Algorithm} (TA), a {\it convex-hull membership} algorithm that given $b_k=Ax_k$ in the ellipsoid $E_{A,\rho}=\{Ax: \Vert x \Vert \leq \rho\}$, it either computes an improved approximation $b_{k+1}=Ax_{k+1}$ or proves $b \not \in E_{A,\rho}$. We then give a dynamic variant of TA, the {\it Centering Triangle Algorithm} (CTA), generating residual, $r_k=b -Ax_k$ via the iteration of $F_1(r)=r-(r^THr/r^TH^2r)Hr$, where $H=AA^T$. If $A$ is symmetric PSD, $H$ can be taken as $A$. Next, for each $t=1, \dots, m$, we derive $F_t(r)=r- \sum_{i=1}^t \alpha_{t,i}(r) H^i r$ whose iterations correspond to a Krylov subspace method with restart. If $\kappa^+(H)$ is the ratio of the largest to smallest positive eigenvalues of $H$, when $Ax=b$ is consistent, in $k=O({\kappa^+(H)}{t^{-1}} \ln \varepsilon^{-1})$ iterations of $F_t$, $\Vert r_k \Vert \leq \varepsilon$. Each iteration takes $O(tN+t^3)$ operations, $N$ the number of nonzero entries in $A$. By directly applying $F_t$ to the normal equation, we get $\Vert A^TAx_k - A^Tb \Vert \leq \varepsilon$ in $O({\kappa^+(AA^T)}{t}^{-1} \ln \varepsilon^{-1})$ iterations. On the other hand, given any residual $r$, we compute $s$, the degree of its minimal polynomial with respect to $H$ in $O(sN+s^3)$ operations. Then $F_s(r)$ gives the minimum-norm solution of $Ax=b$ or an exact solution of $A^TAx=A^Tb$. The proposed algorithms are simple to implementation and theoretically robust. We present sample computational results, comparing the performance of CTA with CG and GMRES methods. The results support CTA as a highly competitive option.
翻译:我们提出了创新算法,用于计算 $Ax=b$ 或{\it 正规方程} $A^TAx=A^Tb$ 的精确或近似(最小范数)解,其中 $A$ 为任意秩的 $m \times n$ 实矩阵。当 $A$ 为对称半正定矩阵时,我们给出了更高效的算法。首先引入{\it 三角形算法}(TA),一种{\it 凸包成员}算法:给定椭球 $E_{A,\rho}=\{Ax: \Vert x \Vert \leq \rho\}$ 中的 $b_k=Ax_k$,该算法要么计算改进近似 $b_{k+1}=Ax_{k+1}$,要么证明 $b \not \in E_{A,\rho}$。随后我们给出TA的动态变体——{\it 中心三角形算法}(CTA),通过迭代 $F_1(r)=r-(r^THr/r^TH^2r)Hr$ 生成残差 $r_k=b -Ax_k$,其中 $H=AA^T$。若 $A$ 对称半正定,则 $H$ 可取为 $A$。进一步,对每个 $t=1, \dots, m$,我们推导出 $F_t(r)=r- \sum_{i=1}^t \alpha_{t,i}(r) H^i r$,其迭代对应带重启的Krylov子空间方法。设 $\kappa^+(H)$ 为 $H$ 的最大与最小正特征值之比,当 $Ax=b$ 一致时,经过 $k=O({\kappa^+(H)}{t^{-1}} \ln \varepsilon^{-1})$ 次 $F_t$ 迭代可得 $\Vert r_k \Vert \leq \varepsilon$。每次迭代需 $O(tN+t^3)$ 次运算,其中 $N$ 为 $A$ 的非零元个数。通过将 $F_t$ 直接应用于正规方程,可在 $O({\kappa^+(AA^T)}{t}^{-1} \ln \varepsilon^{-1})$ 次迭代内实现 $\Vert A^TAx_k - A^Tb \Vert \leq \varepsilon$。另一方面,对任意残差 $r$,我们可在 $O(sN+s^3)$ 次运算内计算其关于 $H$ 的最小多项式次数 $s$,此时 $F_s(r)$ 给出 $Ax=b$ 的最小范数解或 $A^TAx=A^Tb$ 的精确解。所提算法实现简单且理论稳健。我们给出算例结果,将CTA与共轭梯度法和广义最小残差法的性能进行比较。结果支持CTA作为极具竞争力的选项。