Reputation: 33
Sympy functions
I have not used sympy so far, but i feel like trying it out as it seems to be a good to define functions you want to optimize. 
, 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 
Again it's simple when using normal functions, but if anyone can help me code this function in sympy, i would appreciate it! To sum it up, i would like to code this function, and having the flexibility to let the sum run over a vector of the one such as the example i provided. Also, i am basically doing this to optimize this function using something like newton raphson method. It's just a toy example to get a feel of the algorithm. note: theta is an unknown that i will optimize over. Thanks for any hint:)!
Y
Upvotes: 0
Views: 318
Answers (1)
Reputation: 14480
You can certainly use sympy to find the derivative:
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
You can also lambdify these expressions for faster evaluation (note that I made the sum go from 0 to n-1 to match Python indexing):
In [138]: lambdify((theta, x, n), objective)(1, [1,2,3], 3)
Out[138]: -5.736774750542246
You can also try solving analytically for small input samples e.g.:
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 ⎦
For larger inputs an analytic solution won't work (the equation is polynomial so Abel-Ruffini limits this) but you can solve them numerically using sympy if you want:
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
For large input data you will get better speed using the lambdified objective and derivative functions with something like scipy's minimize:
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])
Note that the result is not as accurate as the one computed using nsolve though.
Also the lambdified function in this case is not as efficient as hand-written numpy code because it has not vectorised the sum:
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))))
Upvotes: 1