I've been working with Moore-Penrose pseudoinverse matrices in the same practical/applied engineering system which I've referred to in a previous question to get solutions where matrix inverses cannot be used.
Doing the calculations via this process involves a least-squares method which is 'built in' to the pseudoinverse matrix. I've learned that the final answer outputted is the solution with the "minimal Euclidean norm" which I believe is the solution which satisfies the initial conditions with the highest associated R-squared value.
I have manually done some calculations and found that there exist solutions that are not the minimal Euclidean norm solution but which are practically fine, or in some cases better to use in my experiments. I good analogy here I think is that the minimal Euclidean norm may give an answer as 2.99999999999 but a more practical and contextually better answer for me would be 2.9.
As an example, I have the following matrices:
Which, when using the pseudoinverse or via the RREF/Gaussian method I get the following Matrix X:
But this is one of many Matrix Xs that satisfy the equation? How can I distinguish between solutions? I want to practically choose the best solution.
Note also the system in RREF form which shows that the system indeed has multiple solutions.
My question is: how can I visualise the multiple solutions for my practical system so that I can manually pick the better solutions which makes the most realistic sense?
Here is the source code with a detailed description
#These are the modules used for this guide
from sympy.interactive import printing
# This allows for a neater display using Latex
printing.init_printing(use_latex = True)
from sympy import Eq, solve_linear_system, Matrix
import numpy as np
import sympy as sp
from sympy import *
import param
import math
from sympy.plotting import plot
from sympy import *
import matplotlib.pyplot as plt
from IPython.display import display, HTML
display(HTML("<style>.container { width:95% !important; }</style>"))
#EXAMPLE 1 (Lit)
# Linear Matrix Equation Ax=B
print('Matrix A')
A = Matrix([(0.00895360376159419, 0.00441046995637796, 0.0166432828542564, 0, 0),
(0.00859661406287653, 0, 0, 0, 0),
((3.74133629860424*(10**-6)), 0.00342051671099893, 0, 0.025001806380511, 0),
(0, 0.0664090945967354, 0, 0, 0.111016870123584),
(0, 0, 0, 0, 0)])
display(A)
detA = A.det()
print('Matrix B')
B = Matrix([2, 1, 0.8, 11, 0])
display(B)
#**************************METHOD 1: Solving for X via Matrix Inverse Method ("X_m1"):**************************
print('**************************#METHOD 1: Solving for X via Matrix Inverse Method ("X_m1")**************************')
if detA != 0:
A_inv = A.inv()
else:
A_inv = A.pinv()
A_invdotB = (Matrix([A_inv.dot(B)])).transpose()
print('Matrix X (Matrix Inverse Method)')
X_m1 = ((A_invdotB) / (sum(A_invdotB)) * 100).applyfunc(lambda x: round(x, 2))
display(X_m1)
print('This is the same and expected answer as my parent model/code')
#**************************METHOD 2: Solving for X via Gaussian Elimination Method ("X_gauss"):**************************
print('**************************#METHOD 2: Solving for X via Gaussian Elimination Method ("X_gauss"):**************************')
System = A.col_insert(100, B)
#display(System)
A_rref = System.rref()[0] #Puts Matrix A into Reduced Row Echelon Form (rref) without showing pivot coordinates
display(A_rref)
print('Note that:\n - A single solution is defined with a normal/standard RREF Matrix.\n - A many-solutions system is characterised with a RREF matrix that has a row of all zeroes. \n - A no-solutions system is defined by a zero row RREFF with the exception of a non-zero B matrix element \n')
print('Matrix X (Matrix RREF Method)')
X_m2 = (A_rref.col(-1) / sum(A_rref.col(-1)) * 100).applyfunc(lambda x: round(x, 2))
display(X_m2)
print('This is the same and expected answer as my parent model/code')
