2016-09-18 218 views
3

這裏的情況是:我有一個符號函數lamb它是變量z的元素和變量h的函數元素的函數。這裏是羔羊象徵功能的圖像如何獲得梯度和Hessian | Sympy

enter image description here

現在我想計算這個函數關於變量eta和xi的梯度和Hessian。當然我使用Google搜索,但是我找不到這樣做的直接方法。我發現的是here,但正如我所說的,它並不是這種情況的最佳方法。任何想法? 貝婁,源代碼。謝謝。

from sympy import Symbol, Matrix, Function, simplify 

eta = Symbol('eta') 
xi = Symbol('xi') 

x = Matrix([[xi],[eta]]) 

h = [Function('h_'+str(i+1))(x[0],x[1]) for i in range(3)] 
z = [Symbol('z_'+str(i+1)) for i in range(3)] 

lamb = 0 
for i in range(3): 
    lamb += 1/(2*sigma**2)*(z[i]-h[i])**2 
simplify(lamb) 

回答

7

您可以使用Stelios建議的Pythonic方式,或者使用SymPy最近添加的一些功能:

In [14]: from sympy.tensor.array import derive_by_array 

In [15]: derive_by_array(lamb, (eta, xi)) 
Out[15]: 
[-(z_1 - h_1(xi, eta))*Derivative(h_1(xi, eta), eta)/sigma**2 - (z_2 - h_2(xi, 
eta))*Derivative(h_2(xi, eta), eta)/sigma**2 - (z_3 - h_3(xi, eta))*Derivativ 
e(h_3(xi, eta), eta)/sigma**2, -(z_1 - h_1(xi, eta))*Derivative(h_1(xi, eta), 
xi)/sigma**2 - (z_2 - h_2(xi, eta))*Derivative(h_2(xi, eta), xi)/sigma**2 - (z 
_3 - h_3(xi, eta))*Derivative(h_3(xi, eta), xi)/sigma**2] 

不幸的是,打印機仍然缺少用於N-暗淡陣列,則可以通過將它們轉換爲一個可視化列表(或,可替換地,使用.tomatrix()):

In [16]: list(derive_by_array(lamb, (eta, xi))) 
Out[16]: 
⎡     ∂        ∂       
⎢ (z₁ - h₁(ξ, η))⋅──(h₁(ξ, η)) (z₂ - h₂(ξ, η))⋅──(h₂(ξ, η)) (z₃ - h₃(ξ, η 
⎢     ∂η        ∂η       
⎢- ──────────────────────────── - ──────────────────────────── - ───────────── 
⎢    2        2        
⎣    σ        σ        

    ∂        ∂        ∂   
))⋅──(h₃(ξ, η)) (z₁ - h₁(ξ, η))⋅──(h₁(ξ, η)) (z₂ - h₂(ξ, η))⋅──(h₂(ξ, η)) 
    ∂η        ∂ξ        ∂ξ   
───────────────, - ──────────────────────────── - ──────────────────────────── 
2        2        2    
σ        σ        σ    

        ∂   ⎤ 
    (z₃ - h₃(ξ, η))⋅──(h₃(ξ, η))⎥ 
        ∂ξ   ⎥ 
- ────────────────────────────⎥ 
       2    ⎥ 
       σ    ⎦ 

對於海森,只需重複兩次:

In [18]: list(derive_by_array(derive_by_array(lamb, (eta, xi)), (eta, xi))) 
Out[18]: 
⎡     2        2       
⎢     ∂        ∂       
⎢ (z₁ - h₁(ξ, η))⋅───(h₁(ξ, η)) (z₂ - h₂(ξ, η))⋅───(h₂(ξ, η)) (z₃ - h₃(ξ, 
⎢     2        2       
⎢     ∂η        ∂η       
⎢- ───────────────────────────── - ───────────────────────────── - ─────────── 
⎢     2        2       
⎣    σ        σ        

     2                  
     ∂       2     2     2  
η))⋅───(h₃(ξ, η)) ⎛∂   ⎞ ⎛∂   ⎞ ⎛∂   ⎞  
     2    ⎜──(h₁(ξ, η))⎟ ⎜──(h₂(ξ, η))⎟ ⎜──(h₃(ξ, η))⎟  (z 
    ∂η    ⎝∂η   ⎠ ⎝∂η   ⎠ ⎝∂η   ⎠  
────────────────── + ─────────────── + ─────────────── + ───────────────, - ── 
    2      2     2     2    
    σ      σ     σ     σ    


       2         2       
       ∂         ∂       
₁ - h₁(ξ, η))⋅─────(h₁(ξ, η)) (z₂ - h₂(ξ, η))⋅─────(h₂(ξ, η)) (z₃ - h₃(ξ, 
       ∂ξ ∂η        ∂ξ ∂η       
───────────────────────────── - ─────────────────────────────── - ──────────── 
       2         2        
      σ         σ        


     2                  
     ∂    ∂   ∂    ∂   ∂    
η))⋅─────(h₃(ξ, η)) ──(h₁(ξ, η))⋅──(h₁(ξ, η)) ──(h₂(ξ, η))⋅──(h₂(ξ, η)) 
    ∂ξ ∂η    ∂η   ∂ξ    ∂η   ∂ξ    
─────────────────── + ───────────────────────── + ───────────────────────── + 
    2        2       2    
    σ        σ       σ    


               2        
∂   ∂         ∂        
──(h₃(ξ, η))⋅──(h₃(ξ, η)) (z₁ - h₁(ξ, η))⋅─────(h₁(ξ, η)) (z₂ - h₂(ξ, η)) 
∂η   ∂ξ        ∂ξ ∂η        
─────────────────────────, - ─────────────────────────────── - ─────────────── 
      2        2         
      σ        σ         


    2         2          
    ∂         ∂    ∂   ∂   
⋅─────(h₂(ξ, η)) (z₃ - h₃(ξ, η))⋅─────(h₃(ξ, η)) ──(h₁(ξ, η))⋅──(h₁(ξ, η)) 
∂ξ ∂η        ∂ξ ∂η    ∂η   ∂ξ   
──────────────── - ─────────────────────────────── + ───────────────────────── 
2         2        2   
σ         σ        σ    


                      ∂ 
    ∂   ∂    ∂   ∂    (z₁ - h₁(ξ, η))⋅── 
    ──(h₂(ξ, η))⋅──(h₂(ξ, η)) ──(h₃(ξ, η))⋅──(h₃(ξ, η))      
    ∂η   ∂ξ    ∂η   ∂ξ        ∂ξ 
+ ───────────────────────── + ─────────────────────────, - ────────────────── 
       2       2        2 
       σ       σ        σ 

2        2        2    
           ∂        ∂    
─(h₁(ξ, η)) (z₂ - h₂(ξ, η))⋅───(h₂(ξ, η)) (z₃ - h₃(ξ, η))⋅───(h₃(ξ, η)) 
2        2        2    
           ∂ξ        ∂ξ    
─────────── - ───────────────────────────── - ───────────────────────────── + 
          2        2     
          σ        σ     

                ⎤ 
       2     2     2⎥ 
⎛∂   ⎞ ⎛∂   ⎞ ⎛∂   ⎞ ⎥ 
⎜──(h₁(ξ, η))⎟ ⎜──(h₂(ξ, η))⎟ ⎜──(h₃(ξ, η))⎟ ⎥ 
⎝∂ξ   ⎠ ⎝∂ξ   ⎠ ⎝∂ξ   ⎠ ⎥ 
─────────────── + ─────────────── + ───────────────⎥ 
     2     2     2  ⎥ 
     σ     σ     σ  ⎦ 
4

可以簡單地計算出梯度向量 「手動」(假設變量被排序爲(z1, z2, z3, eta)):

[lamb.diff(x) for x in z+[eta]] 

類似地,對於Hessian矩陣:

[[lamb.diff(x).diff(y) for x in z+[eta]] for y in z+[eta]]