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: enter image description here Which, when using the pseudoinverse or via the RREF/Gaussian method I get the following Matrix X:

enter image description here

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.

enter image description here

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)])

detA = A.det()

print('Matrix B')
B = Matrix([2, 1, 0.8, 11, 0])

#**************************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()
    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))

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)
A_rref = System.rref()[0] #Puts Matrix A into Reduced Row Echelon Form (rref) without showing pivot coordinates
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))

print('This is the same and expected answer as my parent model/code')

  • 1
    $\begingroup$ For multi solutions, any solution plus any element inside the kernel is also a solution. Visualizing this kernel can be difficult. However, if your criteria for the best solution can be something like “closest to point $x^{*}$ w.r.t. Euclidean distance” then simply modify the problem to $A(x-x^{*})=b$ $\endgroup$
    – acat3
    Commented Dec 14, 2022 at 8:25
  • $\begingroup$ I have strangely found that the ker(A) (or nullspace(A)) is equal to the negative * the elements in the 2nd last column of the RREF form (with the exception of the bottom row element being 1). I'm not sure what the implication of that here is but it just so happens that you also mention the kernel here... could you please elaborate on what you're suggesting so that I can try and input something in my code? $\endgroup$
    – Hendrix13
    Commented Dec 14, 2022 at 9:33
  • $\begingroup$ @RezhaAdrianTanuharja I'd love a more detailed answer here with a worked example if possible, I was never able to work through this. Even a way to quantify the number of solutions would be great (obviously provided they're not infinite). $\endgroup$
    – Hendrix13
    Commented Nov 12, 2023 at 4:56


You must log in to answer this question.