numdifftools: Berechnen einer Jakobi-Matrix

[SimPy]

Hallo Forum,

habe hier die "numdifftools" Bibliothek heruntergeladen und versuche damit die Jakobimatrix einer Funktion zu berechnen. Hier habe ich mal folgenden Miniquellcode:

Code: Alles auswählen

import numpy as numpy
import numdifftools as nd

def fun(x):
  r = x[0]*2 + x[1]*5 + x[2]**2
  return r

Jac = nd.Jacobian(fun) 
print(Jac([1,1,1]))

Eigentlich erwarte ich das allgemeine Ergebnis:
Jac = [[dr/x[0] dr/x[1] dr/x[2]]]
d.h. die Ableitung der Funktion nach jeder einzelnen Variablen.
Meinem Verständnis nach müsste das Ergebnis Folgendes sein:
Jac = [[2 5 2*x[2]]]
Wenn ich nun den Lösungvektor [1,1,1] einsetze in der Form Jac([1,1,1]) hätte ich daher eigentlich als Ergebnis
[[2 5 2]] erwartet. Jedoch erscheint als Ergebnis
[[0.46791994 1.169799986 0.90064399]]
Könnt Ihr mir das erklären?

Und eine andere Frage noch hinterher: Benutzt ihr vielleicht eine andere Methode, die Jakobimatrix zu berechnen? Ich habe keine Standardbibliothek zur automatischen Differenzierung gefunden. Die oben verwendete Bibliothek ist eine finitie Elementemethode, bei der die Lösung lediglich angenähert wird, jedoch nicht exakt berechnet wird. Eine exakte Berechnung wäre natürlich besser.

Danke im Voraus!

[SimPy]

Hallo Python-Forum,

ich habe jetzt die Bibiothek algoPy entdeckt.

Code: Alles auswählen

from algopy import *

def fun(x):
  r = x[0]*2 + x[1]*5 + x[2]**2
  return r

# AlgoPy
x = UTPM.init_jacobian([1,1,1])               # definiere den Vektor, der in die Jakobi-Matrix eingesetzt werden soll
y = fun(x)                                    
algopy_jacobian = UTPM.extract_jacobian(y)
print 'jacobian = ',algopy_jacobian

Und jetzt kommt auch der richtige Vektor [[2 5 2]] heraus. Ich weiß nicht, was ich zuvor falsch gemacht habe, aber ich kann auch mit algoPy sehr gut leben. Ich werde jetzt damit weiterarbeiten, kann jedoch gut sein, dass ich dann noch auf weitere Probleme treffe, die ich mit Euch teile

[SimPy]

Könnt Ihr mir auf die Sprünge helfen, was ich hier falsch mache?

Code: Alles auswählen

import numpy, algopy
from algopy import UTPM, exp

def fun(x):
  return [2*x[0] + 4*x[1] + 6*x[2] + exp(7*x[2] - 4*x[1]), 3*x[0] + 6*x[1] + 9*x[2] + exp(2*x[2] - 5*x[1])]

# ----------------------------------
# AlgoPy: reverse mode using a computational graph
# STEP 1: trace the function evaluation
cg = algopy.CGraph()
x = algopy.Function(numpy.zeros(3))
# x = algopy.Function([1,2,3])
y = fun(x)
cg.trace_off()
cg.independentFunctionList = [x]
cg.dependentFunctionList = [y]

# STEP 2: use the computational graph to evaluate derivatives
print 'Jacobian =', cg.jacobian([1,1,1])
# ----------------------------------

Es erscheint die Fehlermeldung:

Code: Alles auswählen

List object has no attribute 'size'

Verwende ich hingegen

Code: Alles auswählen

def fun(x):
  return numpy.array([2*x[0] + 4*x[1] + 6*x[2] + exp(7*x[2] - 4*x[1]), 3*x[0] + 6*x[1] + 9*x[2] + exp(2*x[2] - 5*x[1])])

erscheint die Fehlermeldung:

Code: Alles auswählen

AttributeError: numpy.ndarray has no attribute 'xbar' 

[BlackJack]

@SimPy: Ein kompletter Traceback wäre hilfreich. So weiss man ja gar nicht *wo* das Problem auftritt.

[SimPy]

Code: Alles auswählen

print 'Jacobian = ', cg.jacobian([1,1,1])
File "C:\Python26\lib\site-packages\algopy\tracer\tracer.py", line 299, in jacobian
M = self.dependentFunctionList[0].size
AttributeError: 'list' object has no attribute 'size'

[SimPy]

@ BlackJack: Ich habe nun eine Lösung gefunden. Die sieht so aus:

Code: Alles auswählen

import numpy, algopy
from algopy import UTPM, exp

def eval_f(x):
    """ some function """
    y = algopy.zeros(2,dtype = x)
    y[0] = x[0]*x[1]*x[2] + exp(x[0])*x[1]
    y[1] = x[0] + x[1] + exp(x[2])
    return y

# forward mode without building the computational graph
# -----------------------------------------------------
x = UTPM.init_jacobian([1,1,1])
y = eval_f(x)
algopy_jacobian = UTPM.extract_jacobian(y)
print 'jacobian = '
print algopy_jacobian

So funktioniert es dann auch. Leider habe ich Probleme mit dem fancy indexing, dass ich gerade in Kombination mit algopy benutzen will. Aber ich werde, weil es mittlerweile nichts mehr mit "numdifftools" zu tun hat, einen neuen Beitrag schreiben.