I'm thinking how to program the evaluation of the expression $A (B+A)^{-1}$, where $A$ is a diagonal matrix with elements equal or greater than zero, and B is a square symmetric definite positive matrix. The difficulty resides in the fact that the elements in the diagonal matrix $A$ can be infinity. If all the elements in $A$ are infinity, then $A (B+A)^{-1} = I$
A straightforward implementation is easy to do, because a standard Cholesky decomposition can be applied for obtaining that inverse, and if you assign high values such as $10^{10}$ instead of infinity, you get the numerical results you would expect. However, I don't like using this technique, because high values instead of infinity can cause accuracy problems, and it's not a tidy solution at all.
I did it symbolically in Mathematica, using 6x6 matrices to begin with (I foresee using it in matrices of up to 40x40 or so), and found that the result of the expression is always a diagonal matrix (I don't have a proof for that, but Mathematica gave me a diagonal matrix for all posible cases in 6x6).
Would I get the desired solution if I rearrange the rows and columns in $A$ and $B$ so that I solve it only for the $n \times n$ submatrix with non-infinity values in $A$, and then take the rest of the result as a unit matrix?
Or do you foresee another way of implementing it?
Note: You can try this example in Mathematica, in which B is a 6x6 symmetric matrix and with only one non-infinite value in A:
b = {{b11, b12, b13, b14, b15, b16}, {b12, b22,
b23, b24, b25, b26}, {b13, b23, b33, b34, b35,
b36}, {b14, b24, b34, b44, b45, b46}, {b15, b25,
b35, b45, b55, b56}, {b16, b26, b36, b46, b56,
b66}}
a = DiagonalMatrix[{v, s, s, s, s, s}]
Limit[a*Inverse[b + a], s -> Infinity, Direction -> 1]
And the solution it gives is:
{{v/(v + b11), 0, 0, 0, 0, 0}, {0, 1, 0, 0, 0, 0}, {0, 0, 1, 0, 0,
0}, {0, 0, 0, 1, 0, 0}, {0, 0, 0, 0, 1, 0}, {0, 0, 0, 0, 0, 1}}