Accurate numerical solution for structured M-matrix algebraic Riccati equations

This paper is concerned with a M-matrix algebraic Riccati equation (MARE) XDX−AX−XB+C=0 for which A is block-diagonal and its defining matrix W=B−D−CA is a nonsingular or irreducible singular M-matrix. Such an MARE can be decomposed into many coupled algebraic Riccati equations (AREs) that can be solved by the Jacobi- or Gauss–Seidel-like iteration updating scheme at the outer-loop while by a doubling algorithm in the inner loop for each coupled ARE, as first proposed by Meini (2013). The goals of this paper are two-fold. One is to resolve a critical technical detail in Meini’s algorithm that was not addressed. It is about whether each ARE in the inner loop has a minimal nonnegative solution. It is proved that the defining matrix of each coupled ARE during a doubling iteration is indeed a nonsingular or irreducible singular M-matrix and, as a result, they do have minimal nonnegative solutions and a doubling algorithm is an efficient way to compute them. The other goal is to design a highly accurate implementation of the doubling algorithm for the inner loop so that all entries of the minimal nonnegative solution to the original MARE are calculated with high entrywise relative accuracies, regardless of their magnitudes. This is made possible by a novel way of constructing triplet representations for the coupled ARE during doubling iterations. Numerical examples are presented to demonstrate that the resulting algorithm can indeed deliver an entrywise relatively accurate solution.