How to obtain a convergent solution iteratively for a linear system of equations?

Question

I am working on a problem that requires an iterative procedure to solve a linear system of equations, the system of equations in matrix form is:

$$\underbrace{\begin{bmatrix} r_{11} & r_{12} & r_{13} & \cdots & r_{1j} \\ r_{21} & r_{22} & r_{23} & \cdots & r_{2j} \\ r_{31} & r_{32} & r_{33} & \cdots & r_{3j} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ r_{i1} & r_{i2} & r_{i3} & \cdots & r_{ij} \end{bmatrix}}_\textit{R} \underbrace{\begin{bmatrix} a_{1} & 0 & 0 & \cdots & 0 \\ 0 & a_{2} & 0 & \cdots & 0 \\ 0 & 0 & a_{3} & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & a_{j} \end{bmatrix}}_\textit{A} + \underbrace{\begin{bmatrix} b_{1} & b_{2} & b_{3} & \cdots & b_{j} \\ b_{1} & b_{2} & b_{3} & \cdots & b_{j} \\ b_{1} & b_{2} & b_{3} & \cdots & b_{j} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ b_{1} & b_{2} & b_{3} & \cdots & b_{j} \end{bmatrix}}_\textit{B} = \underbrace{\begin{bmatrix} c_{1} & c_{1} & c_{1} & \cdots & c_{1} \\ c_{2} & c_{2} & c_{2} & \cdots & c_{2} \\ c_{3} & c_{3} & c_{3} & \cdots & c_{3} \\ \vdots & \vdots & \vdots & \ddots & \vdots \\ c_{i} & c_{i} & c_{i} & \cdots & c_{i} \end{bmatrix}}_\textit{C}\\ R_{i, j}A_{j, j} + B_{i, j} = C_{i, j} $$

Now matrix R is fully known (input), matrices A, B, C are unknown. I am working on an iterative procedure in which I provide initial guesses for A and B, then calculate C. The iteration is to be carried out till I get a converged value for C (output I require). The code is being developed in Python. So far no luck in choosing the initial guess too.

UPDATE 1

To illustrate the problem I have given an example for a $2 \times 2$ matrix:

$$\begin{bmatrix} r_{11} & r_{12} \\ r_{21} & r_{22} \\ \end{bmatrix} \begin{bmatrix} a_{1} & 0 \\ 0 & a_{2} \\ \end{bmatrix} + \begin{bmatrix} b_{1} & b_{2} \\ b_{1} & b_{2} \\ \end{bmatrix} = \begin{bmatrix} c_{1} & c_{1} \\ c_{2} & c_{2} \\ \end{bmatrix}\\$$ this gives us:

$$a_{1}r_{11} + b_{1} = a_{2}r_{12} + b_{2} = c_{1}\\ a_{1}r_{21} + b_{1} = a_{2}r_{22} + b_{2} = c_{2}$$

From these equations I get:

$$a_{2} = a_{1}\frac{(r_{11} - r_{21})}{(r_{12} - r_{22})}; \quad b_{2}=b_{1} + a_{1}\frac{(r_{11} - r_{21})(r_{11} - r_{12})}{(r_{12} - r_{22})}$$

This helps in the iterative method, but I would like to generalize it for a bigger matrix.

UPDATE 2

I have added two reproducible Python codes that should be helpful for a $2 \times 2$ matrix:

Code 1

import numpy as np
R = np.matrix([[2.5, 2.9], [2.3, 2.7]])
m = R.shape[0]
n = R.shape[1]
A = np.matrix(np.zeros(shape=(n,n)))
B = np.matrix(np.zeros(shape=(m,n)))
A[0,0] = 1
B[0,0] = 1
k = 0
for i in range (1, 10000):
    A[0,0] = A[0, 0] - 0.0000001
    B[0,0] = B[0, 0] - 0.0000001
    for j in range(1, n):
        A[j, j] = A[(j-1), (j-1)] * ((R[0, 0] - R[1, 0]) / (R[0, 1] - R[1, 1]))
        B[0, j] = B[0, 0] + A[0, 0] * (((R[0, 0] - R[1, 0]) / (R[0, 1] - R[1, 1])) * (R[0, 0] - R[0, 1]))
        B[1, 0] = B[0, 0]
        B[1, j] = B[0, 0] + A[0, 0] * (((R[0, 0] - R[1, 0]) / (R[0, 1] - R[1, 1])) * (R[0, 0] - R[0, 1]))
    C = R * A + B
    C_convergence = np.all(np.all(np.diff(C, axis = 1) == 0, axis = 1), axis = 0)
    k = k + 1
    print k
print C

Code 2

import numpy as np
R = np.matrix([[2.5, 2.9], [2.3, 2.7]])
m = R.shape[0]
n = R.shape[1]
A = np.matrix(np.zeros(shape=(n,n)))
B = np.matrix(np.zeros(shape=(m,n)))
A[0,0] = 1
B[0,0] = 1
k = 0
C_convergence = False
k = 0
while C_convergence == False:
    A[0,0] = A[0, 0] - 0.0000001
    B[0,0] = B[0, 0] - 0.0000001
    for j in range(1, n):
        A[j, j] = A[(j-1), (j-1)] * ((R[0, 0] - R[1, 0]) / (R[0, 1] - R[1, 1]))
        B[0, j] = B[0, 0] + A[0, 0] * (((R[0, 0] - R[1, 0]) / (R[0, 1] - R[1, 1])) * (R[0, 0] - R[0, 1]))
        B[1, 0] = B[0, 0]
        B[1, j] = B[0, 0] + A[0, 0] * (((R[0, 0] - R[1, 0]) / (R[0, 1] - R[1, 1])) * (R[0, 0] - R[0, 1]))
    C = R * A + B
    C_convergence = np.all(np.all(np.diff(C, axis = 1) == 0, axis = 1), axis = 0)
    k = k + 1
    print k
print C

UPDATE 3

I have written the matrices in an alternative way, which is of the form $AX = B$:

$$\begin{bmatrix} r_{11} & 0 & 0 & \cdots & 0 & 1 & 0 & 0 & \cdots & 0 & -1 & 0 & 0 & \cdots & 0 \\ 0 & r_{12} & 0 & \cdots & 0 & 0 & 1 & 0 & \cdots & 0 & -1 & 0 & 0 & \cdots & 0 \\ 0 & 0 & r_{13} & \cdots & 0 & 0 & 0 & 1 & \cdots & 0 & -1 & 0 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & r_{1j} & 0 & 0 & 0 & \cdots & 1 & -1 & 0 & 0 & \cdots & 0 \\ r_{21} & 0 & 0 & \cdots & 0 & 1 & 0 & 0 & \cdots & 0 & 0 & -1 & 0 & \cdots & 0 \\ 0 & r_{22} & 0 & \cdots & 0 & 0 & 1 & 0 & \cdots & 0 & 0 & -1 & 0 & \cdots & 0 \\ 0 & 0 & r_{23} & \cdots & 0 & 0 & 0 & 1 & \cdots & 0 & 0 & -1 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & r_{2j} & 0 & 0 & 0 & \cdots & 1 & 0 & -1 & 0 & \cdots & 0 \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots &\\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots &\\ r_{i1} & 0 & 0 & \cdots & 0 & 1 & 0 & 0 & \cdots & 0 & 0 & 0 & 0 & \cdots & -1 \\ 0 & r_{i2} & 0 & \cdots & 0 & 0 & 1 & 0 & \cdots & 0 & 0 & 0 & 0 & \cdots & -1 \\ 0 & 0 & r_{i3} & \cdots & 0 & 0 & 0 & 1 & \cdots & 0 & 0 & 0 & 0 & \cdots & -1 \\ \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \ddots & \vdots \\ 0 & 0 & 0 & \cdots & r_{ij} & 0 & 0 & 0 & \cdots & 1 & 0 & 0 & 0 & \cdots & -1 \\ \end{bmatrix} \begin{bmatrix} a_{1} \\ a_{2} \\ a_{3} \\ \vdots \\ a_{j} \\ b_{1} \\ b_{2} \\ b_{3} \\ \vdots \\ b_{j} \\ c_{1} \\ c_{2} \\ c_{3} \\ \vdots \\ c_{i} \\ \end{bmatrix} = 0$$

Now how do I solve this using an iterative method?

Answer 1

So the problem is you want to satisfy $$ R \operatorname{diag}(a) = c\mathbf{1}^T - \mathbf{1}b^T $$ for a given matrix $R$, and unknown column vectors $a,b,c$, and $\mathbf{1}$ is the vector of all $1$'s. This says that $R$ is of rank $2$.

If you know for a fact that $R$ has this form already, then you can compute its singular value decomposition $R = U\Sigma V^T$, or in dyadic form, $$ R = \sum_{i} \sigma_i u_iv_i^T $$ The singular values will be zero beyond the first two, so $$ R = \sigma_1 u_1 v_1^T + \sigma_2 u_2 v_2^T = [\operatorname{diag}(a)^{-1} c]\mathbf{1}^T - [\operatorname{diag}(a) \mathbf{1} ]b^T $$ This allows you to identify all the unknowns directly; one of your right singular vectors should all have the same elements, and the other right singular vector will be $b$. Then $a$ is the left singular vector corresponding to $b$ times minus the singular value. With the left singular vector corresponding to the constant right singular vector, and the obtained $a$, you can back out $c$.