Created
December 14, 2018 21:44
-
-
Save Upabjojr/290b720f912dd6c15b87e2e8a47ce693 to your computer and use it in GitHub Desktop.
Script to check results of derivatives of matricial expression
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
from symengine import * | |
X = Matrix(var("X1:10")).reshape(3, 3) | |
A = Matrix(var("A1:10")).reshape(3, 3) | |
B = Matrix(var("B1:10")).reshape(3, 3) | |
C = Matrix(var("C1:10")).reshape(3, 3) | |
D = Matrix(var("D1:10")).reshape(3, 3) | |
a = Matrix(var("a1:4")).reshape(3, 1) | |
b = Matrix(var("b1:4")).reshape(3, 1) | |
c = Matrix(var("c1:4")).reshape(3, 1) | |
def Trace(x): | |
return x[0,0] + x[1,1] + x[2,2] | |
def Inverse(x): | |
return Matrix([[(X[1, 1]*X[2, 2] - X[1, 2]*X[2, 1])/(X[0, 0]*X[1, 1]*X[2, 2] - X[0, 0]*X[1, 2]*X[2, 1] - X[0, 1]*X[1, 0]*X[2, 2] + X[0, 1]*X[1, 2]*X[2, 0] + X[0, 2]*X[1, 0]*X[2, 1] - X[0, 2]*X[1, 1]*X[2, 0]), (-X[0, 1]*X[2, 2] + X[0, 2]*X[2, 1])/(X[0, 0]*X[1, 1]*X[2, 2] - X[0, 0]*X[1, 2]*X[2, 1] - X[0, 1]*X[1, 0]*X[2, 2] + X[0, 1]*X[1, 2]*X[2, 0] + X[0, 2]*X[1, 0]*X[2, 1] - X[0, 2]*X[1, 1]*X[2, 0]), (X[0, 1]*X[1, 2] - X[0, 2]*X[1, 1])/(X[0, 0]*X[1, 1]*X[2, 2] - X[0, 0]*X[1, 2]*X[2, 1] - X[0, 1]*X[1, 0]*X[2, 2] + X[0, 1]*X[1, 2]*X[2, 0] + X[0, 2]*X[1, 0]*X[2, 1] - X[0, 2]*X[1, 1]*X[2, 0])], [(-X[1, 0]*X[2, 2] + X[1, 2]*X[2, 0])/(X[0, 0]*X[1, 1]*X[2, 2] - X[0, 0]*X[1, 2]*X[2, 1] - X[0, 1]*X[1, 0]*X[2, 2] + X[0, 1]*X[1, 2]*X[2, 0] + X[0, 2]*X[1, 0]*X[2, 1] - X[0, 2]*X[1, 1]*X[2, 0]), (X[0, 0]*X[2, 2] - X[0, 2]*X[2, 0])*X[0, 0]/((X[0, 0]*X[1, 1] - X[0, 1]*X[1, 0])*(X[0, 0]*X[2, 2] - X[0, 2]*X[2, 0]) - (X[0, 0]*X[1, 2] - X[0, 2]*X[1, 0])*(X[0, 0]*X[2, 1] - X[0, 1]*X[2, 0])), -(X[0, 0]*X[1, 2] - X[0, 2]*X[1, 0])*X[0, 0]/((X[0, 0]*X[1, 1] - X[0, 1]*X[1, 0])*(X[0, 0]*X[2, 2] - X[0, 2]*X[2, 0]) - (X[0, 0]*X[1, 2] - X[0, 2]*X[1, 0])*(X[0, 0]*X[2, 1] - X[0, 1]*X[2, 0]))], [(X[1, 0]*X[2, 1] - X[1, 1]*X[2, 0])/(X[0, 0]*X[1, 1]*X[2, 2] - X[0, 0]*X[1, 2]*X[2, 1] - X[0, 1]*X[1, 0]*X[2, 2] + X[0, 1]*X[1, 2]*X[2, 0] + X[0, 2]*X[1, 0]*X[2, 1] - X[0, 2]*X[1, 1]*X[2, 0]), -(X[0, 0]*X[2, 1] - X[0, 1]*X[2, 0])*X[0, 0]/((X[0, 0]*X[1, 1] - X[0, 1]*X[1, 0])*(X[0, 0]*X[2, 2] - X[0, 2]*X[2, 0]) - (X[0, 0]*X[1, 2] - X[0, 2]*X[1, 0])*(X[0, 0]*X[2, 1] - X[0, 1]*X[2, 0])), (X[0, 0]*X[1, 1] - X[0, 1]*X[1, 0])*X[0, 0]/((X[0, 0]*X[1, 1] - X[0, 1]*X[1, 0])*(X[0, 0]*X[2, 2] - X[0, 2]*X[2, 0]) - (X[0, 0]*X[1, 2] - X[0, 2]*X[1, 0])*(X[0, 0]*X[2, 1] - X[0, 1]*X[2, 0]))]]) | |
expr = Trace(Inverse(X.T*C*X)*A) | |
part1 = -(C*X*Inverse(X.T*C*X)) | |
part2 = Inverse(X.T*C*X) | |
result = part1*A*part2 + part1*A.T*part2 | |
assert expr.diff(X) == part1*A*part2 + part1*A.T*part2 | |
def verify_result(expr, result): | |
for i in range(3): | |
for j in range(3): | |
v = X[i, j] | |
assert expr.diff(v) == result[i, j] |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
@asmeurer this script can help check matrix derivatives by explicitly replacing them with 3x3 matrices. Unfortunately it's extremely slow.