同情函数

Question

到目前为止，我还没有使用过 sympy，但我想尝试一下，因为定义要优化的函数似乎很好。 , this is the function that i am trying to code with sympy. I know i could do it simply using a normal function, but it would be cool for example to code it in sympy and then it would be quite easy i would think for python to find the derivative for example. I want the sum to run over 同样，使用普通函数时它很简单，但是如果有人可以帮助我在 sympy 中编写此函数，我将不胜感激！总而言之，我想编写此函数的代码，并可以灵活地让总和运行超过我提供的示例的向量。另外，我这样做基本上是为了使用类似 newton raphson 方法的方法来优化此功能。这只是一个用来感受算法的玩具示例。注意：theta 是一个未知数，我将对其进行优化。感谢您的任何提示:)！

是

Answer 1

当然可以用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, 0, 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

您还可以对这些表达式进行 lambdify 以加快计算速度（请注意，我将总和从 0 变为 n-1 以匹配 Python 索引）：

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

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

In [142]: vec = [1, 2, 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, 2, 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

对于大型输入数据，使用 lambdified objective 和衍生函数以及类似 scipy 的最小化：

可以获得更快的速度

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

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

In [160]: from scipy.optimize import minimize                                                                                                  

In [161]: minimize(f, 1, jac=fp, args=([1, 2, 3], 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 函数不如 hand-written 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))))

同情函数

Sympy functions

sympy

mathematical-optimization