Sympy功能

问题描述

到目前为止，我还没有使用过sympy，但是我想尝试一下，因为定义要优化的函数似乎很好。

，这是我尝试使用sympy进行编码的功能。我知道我可以简单地使用常规函数来完成此操作，但是例如以sympy对其进行编码会很酷，然后我会很容易地想到python查找派生类。我希望总和超过

同样，使用常规功能很简单，但是如果有人可以帮助我用sympy编写此功能，我将不胜感激！总结起来，我想对这个函数进行编码，并具有让总和在一个向量上运行的灵活性，例如我提供的示例。另外，我基本上是使用牛顿拉夫森（Newton Raphson）方法来优化此功能的。这只是一个玩具示例，可以让您对该算法有所了解。注意：theta是一个未知数，我将对其进行优化。感谢您提供任何提示：）！

解决方法

您当然可以使用sympy来找到导数：

In [132]: x = IndexedBase('x',real=True)                                                                                                      

In [133]: n,i = symbols('n,i',integer=True,positive=True)                                                                                  

In [134]: theta = Symbol('theta',real=True)                                                                                                   

In [135]: objective = -n*log(pi) - Sum(log(1 + (x[i] - theta)**2),(i,n-1))                                                                

In [136]: objective                                                                                                                            
Out[136]: 
            n - 1                      
             ___                       
             ╲                         
              ╲      ⎛           2    ⎞
-n⋅log(π) -   ╱   log⎝(-θ + x[i])  + 1⎠
             ╱                         
             ‾‾‾                       
            i = 0                      

In [137]: objective.diff(theta)                                                                                                                
Out[137]: 
 n - 1                 
  ____                 
  ╲                    
   ╲     2⋅θ - 2⋅x[i]  
    ╲  ────────────────
-   ╱             2    
   ╱   (-θ + x[i])  + 1
  ╱                    
  ‾‾‾‾                 
 i = 0

您还可以对这些表达式进行lambd化，以便进行更快的评估（请注意，我将总和从0改为n-1以匹配Python索引）：

In [138]: lambdify((theta,x,n),objective)(1,[1,2,3],3)                                                                                    
Out[138]: -5.736774750542246

您也可以尝试分析较小的输入样本，例如：

In [142]: vec = [1,3]                                                                                                                      

In [143]: objective.diff(theta).subs(n,3).doit().subs({x[i]:vec[i] for i in range(3)})                                                        
Out[143]: 
    2⋅θ - 6        2⋅θ - 4        2⋅θ - 2   
- ──────────── - ──────────── - ────────────
         2              2              2    
  (3 - θ)  + 1   (2 - θ)  + 1   (1 - θ)  + 1

In [144]: solve(_,theta)                                                                                                                      
Out[144]: 
⎡           _____________          _____________          _____________          _____________⎤
⎢          ╱   1   √11⋅ⅈ          ╱   1   √11⋅ⅈ          ╱   1   √11⋅ⅈ          ╱   1   √11⋅ⅈ ⎥
⎢2,2 -   ╱  - ─ - ─────,2 +   ╱  - ─ - ─────,2 -   ╱  - ─ + ─────,2 +   ╱  - ─ + ───── ⎥
⎣       ╲╱     3     3         ╲╱     3     3         ╲╱     3     3         ╲╱     3     3   ⎦

对于较大的输入，解析解决方案将不起作用（该方程式是多项式，因此Abel-Ruffini对此进行了限制），但您可以根据需要使用sympy对其进行数值求解：

In [150]: vec = [1,3,4,5]                                                                                                                

In [151]: objective.diff(theta).subs(n,5).doit().subs({x[i]:vec[i] for i in range(5)})                                                        
Out[151]: 
    2⋅θ - 10       2⋅θ - 8        2⋅θ - 6        2⋅θ - 4        2⋅θ - 2   
- ──────────── - ──────────── - ──────────── - ──────────── - ────────────
         2              2              2              2              2    
  (5 - θ)  + 1   (4 - θ)  + 1   (3 - θ)  + 1   (2 - θ)  + 1   (1 - θ)  + 1

In [152]: nsolve(_,theta,1)                                                                                                                  
Out[152]: 3.00000000000000

对于大型输入数据，使用经过层层化的目标函数和派生函数（如scipy's minimum）可以得到更快的速度：

In [158]: f = lambdify((theta,-objective)                                                                                              

In [159]: fp = lambdify((theta,-objective.diff(theta))                                                                                 

In [160]: from scipy.optimize import minimize                                                                                                  

In [161]: minimize(f,1,jac=fp,args=([1,3))                                                                                          
Out[161]: 
      fun: 4.820484018668093
 hess_inv: array([[0.50002255]])
      jac: array([9.77560778e-08])
  message: 'Optimization terminated successfully.'
     nfev: 4
      nit: 3
     njev: 4
   status: 0
  success: True
        x: array([2.00000005])

请注意，结果并不如使用nsolve计算的结果准确。

在这种情况下，lambdified函数的效率不如手写的numpy代码，因为它没有向量化和：

In [169]: print(inspect.getsource(f))                                                                                                          
def _lambdifygenerated(theta,Dummy_167,n):
    return (n*log(pi) + (builtins.sum(log((-theta + Dummy_167[i])**2 + 1) for i in range(0,n - 1+1))))

mathematical-optimization sympy