3
$\begingroup$

I am trying to obtain the solution to a modified linear system. I am comparing two methods to solve this modified linear system, and I'm noticing some issues with one of the methods.

A linear system is written as $$ Ax=b $$

where $A \in \mathbb{R}^{n\times n}$ and $x,b \in \mathbb{R}^n$. For a non-singular $A$ and a known $b$, we can obtain $x$ using $$ x=A^{-1}b $$

I came across a problem where I am basically given something that looks like $Ax=b$, but 1 element in $x$ is known beforehand, and 1 element in $b$ is unknown beforehand.

Suppose the known element in $x$ is $x_i$ and unknown element in $b$ is $b_j$. $b_p$ is known for $ p \in \{1,\ldots,n\}\setminus \{j\}$. $x_q$ is unknown for $q \in \{1,\ldots,n\}\setminus \{i\}$

I want to formulate the system $A'x'=b'$ where $x'$ is the vector of unknowns and $b'$ is the vector of knowns. I want to solve for $x'$.

I've basically identified two ways to solve to do this. Method 1, which does some "swapping" PRIOR to taking an inverse, and Method 2, which does the "swapping" AFTER taking an inverse.

Method 1:

We can write the original system as

$$ \begin{bmatrix} A_{1,1}x_1 + \cdots + A_{1,i-1}x_{i-1} + A_{1,i}x_i + A_{1,i+1}x_{i+1} + \cdots + A_{1,n}x_{n} \\ \vdots\\ A_{j-1,1}x_1 + \cdots + A_{j-1,i-1}x_{i-1} + A_{j-1,i}x_i + A_{j-1,i+1}x_{i+1} + \cdots + A_{j-1,n}x_{n} \\ A_{j,1}x_1 + \cdots + A_{j,i-1}x_{i-1} + A_{j,i}x_i + A_{j,i+1}x_{i+1} + \cdots + A_{j,n}x_{n} \\ A_{j+1,1}x_1 + \cdots + A_{j+1,i-1}x_{i-1} + A_{j+1,i}x_i + A_{j+1,i+1}x_{i+1} + \cdots + A_{j+1,n}x_{n} \\ \vdots\\ A_{n,1}x_1 + \cdots + A_{n,i-1}x_{i-1} + A_{n,i}x_i + A_{n,i+1}x_{i+1} + \cdots + A_{n,n}x_{n} \\ \end{bmatrix} = \begin{bmatrix} b_1 \\ \vdots \\ b_{j-1} \\ b_j \\ b_{j+1} \\ \vdots \\ b_n \end{bmatrix} $$

We can rearrange the terms so that the RHS consists only of knowns and all the known $x$'s are on the RHS (I call this step the "swapping"): $$ \begin{bmatrix} A_{1,1}x_1 + \cdots + A_{1,i-1}x_{i-1} + 0 + A_{1,i+1}x_{i+1} + \cdots + A_{1,n}x_{n} \\ \vdots\\ A_{j-1,1}x_1 + \cdots + A_{j-1,i-1}x_{i-1} + 0 + A_{j-1,i+1}x_{i+1} + \cdots + A_{j-1,n}x_{n} \\ A_{j,1}x_1 + \cdots + A_{j,i-1}x_{i-1} - b_j + A_{j,i+1}x_{i+1} + \cdots + A_{j,n}x_{n} \\ A_{j+1,1}x_1 + \cdots + A_{j+1,i-1}x_{i-1} + 0 + A_{j+1,i+1}x_{i+1} + \cdots + A_{j+1,n}x_{n} \\ \vdots\\ A_{n,1}x_1 + \cdots + A_{n,i-1}x_{i-1} + 0 + A_{n,i+1}x_{i+1} + \cdots + A_{n,n}x_{n} \\ \end{bmatrix} = \begin{bmatrix} b_1 - A_{1,i}x_i \\ \vdots \\ b_{j-1} - A_{j-1,i}x_i \\ -A_{j,i}x_i \\ b_{j+1} - A_{j+1,i}x_i \\ \vdots \\ b_n - A_{n,i}x_i \end{bmatrix} $$

So now we can form $A'x'=b'$ where $$ A' = \begin{bmatrix} A_{1,1} & \cdots & A_{1,i-1} & 0 & A_{1,i+1} & \cdots & A_{1,n} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ A_{j-1,1} & \cdots & A_{j-1,i-1} & 0 & A_{j-1,i+1} & \cdots & A_{j-1,n} \\ A_{j,1} & \cdots & A_{j,i-1} & -1 & A_{j,i+1} & \cdots & A_{j,n} \\ A_{j+1,1} & \cdots & A_{j+1,i-1} & 0 & A_{j+1,i+1} & \cdots & A_{j+1,n} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ A_{n,1} & \cdots & A_{n,i-1} & 0 & A_{n,i+1} & \cdots & A_{n,n} \\ \end{bmatrix} \\ x' = \begin{bmatrix} x_1 \\ \vdots \\ x_{i-1} \\ b_j \\ x_{i+1}\\ \vdots \\ x_n \end{bmatrix} \\ b' = \begin{bmatrix} b_1 - A_{1,i}x_i \\ \vdots \\ b_{j-1} - A_{j-1,i}x_i \\ -A_{j,i}x_i \\ b_{j+1} - A_{j+1,i}x_i \\ \vdots \\ b_n - A_{n,i}x_i \end{bmatrix} $$

So we can now solve for $x'=A'^{-1}b'$. But I am finding that $A'$ has a large condition number in some of my application cases. In one case the condition number was $O(10^{16})$, which is basically numerically singular.

Method 2:

We work with the original $Ax=b$. We take the inverse of $A$, and write $x=A^{-1}b = Cb$ where $C=A^{-1}$. Recall that $b_j$ is unknown and $x_i$ is known.

So if we expand out the terms we can write something that looks like this (note the "swapping" in this method occurs after we have taken the inverse): $$ \begin{bmatrix} x_1 - b_j*C_{1,j}\\ \vdots \\ x_{i-1}-b_j*C_{i-1,j} \\ -b_j*C_{i,j} \\ x_{i+1}-b_j*C_{i+1,j} \\ \vdots\\ x_n-b_j*C_{n,j} \end{bmatrix} = \begin{bmatrix} C_{1,1} & \cdots & C_{1,j-1} & 0 & C_{1,j+1} & \cdots & C_{1,n} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ C_{i-1,1} & \cdots & C_{i-1,j-1} & 0 & C_{i-1,j+1} & \cdots & C_{i-1,n} \\ C_{i,1} & \cdots & C_{i,j-1} & -1 & C_{i,j+1} & \cdots & C_{i,n} \\ C_{i+1,1} & \cdots & C_{i+1,j-1} & 0 & C_{i+1,j+1} & \cdots & C_{i+1,n} \\ \vdots & \vdots & \vdots & \vdots & \vdots & \vdots & \vdots \\ C_{n,1} & \cdots & C_{n,j-1} & 0 & C_{n,j+1} & \cdots & C_{n,n} \\ \end{bmatrix} \begin{bmatrix} b_1 \\ \vdots \\ b_{j-1} \\ x_i \\ b_{j+1} \\ \vdots \\ b_{n} \end{bmatrix} $$

Multiplying $Cb$ gives me the LHS of the above equation. Then I can solve for $b_j$ and use that to get the unknown $x$'s.

Question/Confusion

Both methods should arrive at the same answer (and it does in most cases), but the issue is in Method 1, the matrix that is being inverted, $A'$, is close to being singular in some cases, whereas the matrix $A$ in Method 2 is much more well conditioned. For example, in one case that I conducted, $A'$ had a condition number that was $O(10^{17})$, but $A$'s condition number was only $O(10^6)$. But both methods should arrive at the same answer, so my confusion is why is method 1 so ill conditioned (for some cases) when they're effectively solving the same problem, just with some different sequence of linear algebra operations?

Edit (Including the referenced matrix that is giving me issues)

Here is the matrix that is causing me issues ($n=6$):

$$ A=\begin{bmatrix} 567e+009 & -49e+003 & 0 & 0 & 77e+000 & 623e+006\\ 0 & 0 & 567e+009 & -49e+003 & 260e+000 & -10.1493e+012\\ 590e+003 & 5e+012 & 2e+006 & -77e+015 & 17e+015 & -1e+009\\ -217e+000 & 20e+012 & -77e+000 & -623e+006 & 0 & 0\\ -590e+003 & -4e+012 & -6e+006 & 618e+015 & -2e+015 & -62e+003\\ 0 & 0 & -283e+009 & -8e+000 & -217e+000 & 20e+012 \end{bmatrix} $$ This matrix has a condition number that is $O(10^6)$. This matrix is obviously not that well conditioned to begin with.

For this matrix, $i = 3$ and $j=2$, so we can write $$ $$ A'=\begin{bmatrix} 567e+009 & -49e+003 & 0 & 0 & 77e+000 & 623e+006\\ 0 & 0 & -1 & -49e+003 & 260e+000 & -10.1493e+012\\ 590e+003 & 5e+012 & 0 & -77e+015 & 17e+015 & -1e+009\\ -217e+000 & 20e+012 & 0 & -623e+006 & 0 & 0\\ -590e+003 & -4e+012 & 0 & 618e+015 & -2e+015 & -62e+003\\ 0 & 0 & 0 & -8e+000 & -217e+000 & 20e+012 \end{bmatrix} $$

This matrix has a condition number of about $O(10^{16})$. Matlab is actually telling me that this matrix is rank 5..

MATLAB Script Below (in reference to above matrices):

clc;  clear all; close all
format shorteng

%note this matrix was copied over from a simulation code for debugging
A = [   567322834646.66357       -49079.870115227313        0.0000000000000000        0.0000000000000000        77.820576525741828        623686024.79204798     
        0.0000000000000000        0.0000000000000000        567322834646.66357       -49079.870115227313        260.93116192028197       -10149285653847.125     
        590586.40976551455        4733201765487.0498        1980226.8884948289       -77023718450704160.        17221828521465004.       -1489883811.0520439     
        -217.85980585264346        20382594821474.887       -77.820576525741828       -623686024.79204798        0.0000000000000000        0.0000000000000000     
        -590586.40976551455       -4733201765487.0498       -6613420.0652273158        6.1874039153838144E+017  -2152728565180286.5       -62487.145455156555     
        0.0000000000000000        0.0000000000000000       -283661417322.95770       -8.2338259132967249       -217.85980585264346        20382594821474.887 ]
 
i = 3;
j = 2;

%This is A' on math.stackexchange
A_prime = A;
A_prime(:,i) = 0;
A_prime(j,i) = -1;

fprintf('The rank of A is: %i \n', rank(A))
fprintf('The rank of A_prime is: %i \n', rank(A_prime))
fprintf('The condition number of A is: %e \n', cond(A))
fprintf('The condition number of A_prime is: %e \n', cond(A_prime))    

If you execute this code, you will get the following output:

The rank of A is: 6 
The rank of A_prime is: 5 
The condition number of A is: 1.635215e+06 
The condition number of A_prime is: 8.141279e+17 

So $x'$ is not solvable if using method 1. Or I should say that it appears method 1 is suggesting that $x'$ is not unique.

But if you use method 2, then you can obtain a unique (?) $x'$. All we need to do is invert the original $A$ matrix, do some swapping and some simple algebra to get $x'$.

My gut tells me method 1 is more sound, and that method 2 might be producing a $x'$ vector, but the solution produced is not unique?

$\endgroup$
4
  • $\begingroup$ It seems your matrix for $A'$ does not match what it should based on your equations. The bottom row seems to have a $0$ replacing an element of the matrix that it should not be. Assuming you fix that, does the results change? $\endgroup$
    – spektr
    Mar 7 '19 at 7:36
  • $\begingroup$ I should clarify that I'm talking about your example $A'$ at the bottom of your answer in the Edit section. It doesn't match the $A'$ equation you define earlier, which makes sense. $\endgroup$
    – spektr
    Mar 7 '19 at 14:11
  • $\begingroup$ @spektr I think that’s a typo in the OP. I will check my code when I get to my computer. $\endgroup$
    – 24n8
    Mar 7 '19 at 14:14
  • $\begingroup$ Just checked. It was a typo when I copied it over on here. $\endgroup$
    – 24n8
    Mar 7 '19 at 14:25
1
$\begingroup$

As a toy problem, consider the system of equations

\begin{align} \underbrace{\begin{bmatrix} a & b \\ c & d \end{bmatrix}}_{A} \begin{bmatrix} x_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} r_1 \\ r_2 \end{bmatrix} \end{align}

Let us assume we know $x_1$ and don't know $r_1$. We can rewrite the system of equations to be

\begin{align} \underbrace{\begin{bmatrix} -1 & b \\ 0 & d \end{bmatrix}}_{A'} \begin{bmatrix} r_1 \\ x_2 \end{bmatrix} = \begin{bmatrix} -a x_1 \\ r_2 - c x_1 \end{bmatrix} \end{align}

We know that $\text{det}(A) = ad - bc$. We also readily see that $\text{det}(A') = -d$. Clearly, the determinants between $A$ and $A'$ are very different and we can see that if $|d| \rightarrow 0$, we will begin to have a poorly conditioned system of equations that ultimately become singular. This is essentially Method 1.

If we consider Method 2, we find that

\begin{align} C = A^{-1} &= \begin{bmatrix} a & b \\ c & d\end{bmatrix}^{-1} \\ &= \text{det}(A)^{-1} \begin{bmatrix} d & -b \\ -c & a\end{bmatrix} \end{align}

Using the $C$ found above and separating known and unknown terms in the equation $Cr = x$, we can find the new system of equations to solve is

\begin{align} \underbrace{\text{det}(A)^{-1}\begin{bmatrix} d & -b \\ -c & a \end{bmatrix}}_{C} \begin{bmatrix} r_1 \\ r_2 \end{bmatrix} &= \begin{bmatrix}x_1 \\ x_2 \end{bmatrix} \\ &\Downarrow \\ \underbrace{\begin{bmatrix} \text{det}(A)^{-1} d & 0 \\ -\text{det}(A)^{-1} c & -1 \end{bmatrix}}_{C'} \begin{bmatrix} r_1 \\ x_2 \end{bmatrix} &= \begin{bmatrix}x_1 + \text{det}(A)^{-1} b r_2\\ -\text{det}(A)^{-1}a r_2 \end{bmatrix} \end{align}

We readily know that $\text{det}(C) = \text{det}(A^{-1}) = \text{det}(A)^{-1}$. We can also readily see that $\text{det}(C') = -\text{det}(A)^{-1} d = -\text{det}(C) d$. Interestingly in this case, if $|d| \rightarrow 0$, it is still possible that $\text{det}(C)$ could be large enough to make $\text{det}(C')$ not too small, particularly if $|\text{det}(A)|$ is small. Of course, there's also the chance that $|\text{det}(C)| \ll 1$, which could make some $\text{det}(C') \rightarrow 0$ even if $|d|$ is not too small.

The point here is that Method 2 actually appears to produce a different result if we assume that we compute $C$ as an exact inverse and don't multiply the $Cr = x$ equation by $\text{det}(A)$. Even in arbitrary precision, you could see a difference between using Method 1 versus using Method 2. Taking into account finite precision issues and you may see more dramatic differences.

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.