CODEWITHSUNDEEP
pythonnumpy
suzuichi
suzuichi

Reputation: 13

numpy, same calculation different result

I want to calculate this fomula.
I think the result is A.
So I write a python code using numpy.
But depending on the computation sequence, result is not A.
What brought this on? 

import numpy as np
from numpy import *
from numpy.random import *
import decimal

#generate matrix A
A = randn(180,240)
A = np.array(A, dtype = decimal.Decimal )

#generate matrix P
h,w=A.shape
P = randn(0.9*h,h)
P = np.array(P, dtype = decimal.Decimal )

#it's OK. IA = A
PP = dot(P.T,P)
PPinv = np.linalg.inv(PP)
PPinvPP = dot(PPinv,PP)
PPinvPPinv = np.linalg.inv(PPinvPP)
I = dot(PPinvPPinv,PPinvPP)
IA = dot(I, A)

#I think IA2 must be A. but not A. Why?
PA = dot(P,A)
PPA = dot(P.T,PA)
PPinvPPA = dot(PPinv,PPA)
IA2 = dot(PPinvPPinv, PPinvPPA)

#print result
print "A;%0.2f" % A[0,0]
print "IA:%0.2f" % IA[0,0]
print "IA2:%0.2f" % IA2[0,0]

Upvotes: 0

Views: 206

Answers (3)

Dietrich
Dietrich

Reputation: 5531

Not a direct answer to your question, but you could use Sympy for these kind of problems:

from IPython.display import display
import sympy as sy


sy.init_printing()  # For LaTeX-like pretty printing in IPython

n = 5
A = sy.MatrixSymbol('A', 162, 240)  # dimension 162x240
P = sy.MatrixSymbol('P', 162, 180)  # dimensions 162x180

PTP = P*P.T
ex1 = (PTP.inverse() * PTP).inverse() * PTP.inverse() * PTP * A
display(ex1)  # displays: A

Upvotes: 0

user4999547
user4999547

Reputation:

The main reason is when you use non-square matrix P, where height is less than width, determinant of the PP always has a zero value, but because of a calc error it's != 0. So after this it's impossible to calculate the real PPinv and any forward steps are meaningless.

P = randn(2,3)
P = np.array(P, dtype = decimal.Decimal )

PP = dot(P.T,P)
np.linalg.det(PP) #-5.2536080570332981e-34

So why is IA == A?

I think it's a situation when error*error gives you a normal result.

How to solve it?

  1. Do not use Python for theoretical questions :)
  2. Change P = randn(0.9*h,h) to P = randn(1.1*h,h)

Upvotes: 1

cel
cel

Reputation: 31339

What happens here is quite interesting:

In general your formula is only correct if PP is non-singular.

So why then AI == A?

PP = dot(P.T,P)
PPinv = np.linalg.inv(PP)
PPinvPP = dot(PPinv,PP)
PPinvPPinv = np.linalg.inv(PPinvPP)
I = dot(PPinvPPinv,PPinvPP)
IA = dot(I, A)

There are a couple of things to note here:

  • PP = dot(P.T,P) is singular
  • => PPinv is not a true inverse
  • but PPinvPP is invertible, so I is indeed the identity matrix

Note: You only get AI == A because of your special order of evaluation of the terms.

In the second calculation of IA2 term you don't have this special evaluation order that gives you A as result.

Upvotes: 2

Related Questions