0

I'm trying to run the following code, which is simple matrix algebra:

import numpy as np
from numpy import linalg

A = np.array([[1,2],[1,-1],[1,1]])
b = np.array([[1],[-1],[5]])

r = linalg.inv(A.transpose()@A)@A.transpose()@b

print(r)

The output I obtain is:

[[1.]
[1.]]

However, I'm expecting a strict [[1],[1]], which is what I obtain when doing the calculation by hand. Doing further operations with r in this form starts giving me incorrect results, for example, if I compute A@r I obtain:

[[3.00000000e+00]
 [3.33066907e-16]
 [2.00000000e+00]]

Instead of the expected [[3],[0],[2]]. Is there another way to do matrix algebra to avoid this issue, maybe having python use fractions instead of raw computing for taking the inverse?

shintuku
  • 223
  • 1
  • 6
  • Does this answer your question? [Why does numpy.linalg.solve() offer more precise matrix inversions than numpy.linalg.inv()?](https://stackoverflow.com/questions/31256252/why-does-numpy-linalg-solve-offer-more-precise-matrix-inversions-than-numpy-li) – Kaia Jan 06 '23 at 21:43
  • @Kaia ideally, I'd like to use the inverse operation to confirm whether my hand calculations are correct, instead of solving for an unknown variable. I don't know whether that's possible. That's why I'm not sure the linked question answers my question, since it involves using `np.linalg.solve(A, b)` instead – shintuku Jan 06 '23 at 21:46
  • 3
    If you want to do exact rational arithmetic you'll have to use something like sympy and make a matrix over the rationals. – President James K. Polk Jan 06 '23 at 21:57
  • @PresidentJamesK.Polk I found out just now, sympy gives me exact results. Thank you for the answer! Should I delete this question, or would you like to post the sympy suggestion? it might be of use for someone else wanting linear algebra over rational numbers – shintuku Jan 06 '23 at 21:58
  • 2
    Even better, just answer your own question. – President James K. Polk Jan 06 '23 at 21:59

2 Answers2

2

Turns out that it is better to do matrix algebra with the sympy library for this case. We get our desired (and more readable) result with:

import sympy as sy

A = sy.Matrix([[1,2],[1,-1],[1,1]])
b = sy.Matrix([[1],[-1],[5]])

r = (A.T*A)**-1*A.T*b

print(r)
shintuku
  • 223
  • 1
  • 6
-1

This problem occurs due to floating-point number imprecision. See here and here.

Answer found through looking at this question.

import numpy as np
from numpy import linalg

A = np.array([[1,2],[1,-1],[1,1]])
b = np.array([[1],[-1],[5]])

r = linalg.inv(A.transpose()@A)@A.transpose()@b

print(r.round().astype(int))
print(A@r.round().astype(int))

EDIT: The above answers the question using numpy (which the OP asked for). However, it is better to use decimal but it looks like numpy does not support decimal, but pandas can.

tehCheat
  • 333
  • 2
  • 8
  • 3
    Thanks for the answer! but do we not risk getting inexact results if we actually have fractions in our computations? It seems your approach would round rational numbers – shintuku Jan 06 '23 at 21:51