# [SciPy-user] argmin and optimization

dmitrey openopt@ukr....
Thu Jul 26 07:12:00 CDT 2007

if you use python 2.5 you can use more short:
def P_Omega(y):
from numpy.linalg import norm
n_y = norm(y)
return y if n_y <= 1 else y/n_y
HTH, D.

Nils Wagner wrote:
> dmitrey wrote:
>
>> Nils Wagner wrote:
>>
>>
>>> dmitrey wrote:
>>>
>>>
>>>
>>>> So your Omega is just all x: (x,x)<1?
>>>>
>>>>
>>>>
>>>>
>>> (x,x) <= 1 ("closed ball constraint")
>>>
>>>
>>>
>>>
>>>> (Can you attach non-Latex notation elseware?)
>>>>
>>>>
>>>>
>>>>
>>> I can't see what you mean here.
>>>
>>>
>>>
>> I meant couldn't you describe what does Omega and P_Omega mean w/o Latex?
>> And/or describe howto see Latex strings in Linux?
>>
>>
>>>> Why can't you just set P_Omega(y) = {
>>>> y, for ||y|| <=1
>>>> y/||y||, for ||y|| > 1
>>>> }
>>>> ?
>>>>
>>>>
>>>>
>>>>
>>> I don't see your point.
>>>
>>> Nils
>>>
>>>
>>>
>> I meant is this true that the func P_Omega will yield the result you intend?
>> from numpy.linalg import norm
>> def P_Omega(y):
>>   n_y = norm(y)
>>   if n_y <= 1: return y
>>   else: return y/ny
>> HTH, D
>>
>>
>>>
>>> _______________________________________________
>>> SciPy-user mailing list
>>> SciPy-user@scipy.org
>>> http://projects.scipy.org/mailman/listinfo/scipy-user
>>>
>>>
>>>
>>>
>>>
>>>
>> _______________________________________________
>> SciPy-user mailing list
>> SciPy-user@scipy.org
>> http://projects.scipy.org/mailman/listinfo/scipy-user
>>
>>
> Dmitrey,
>
> Thank you very much for your help. It works fine for me.
> I have attached the script for the sake of completeness.
>
> Cheers,
>                     Nils
>
>
> ------------------------------------------------------------------------
>
> from scipy import *
> from pylab import plot, show
> #
> # Gene H. Golub and Li-Zhi Liao
> # Continous methods for extreme and interior eigenvalue problems
> # Linear Algebra and its Applications Vol. 415 (2006) pp. 31-51
> #
> def merit(x):
>     """ Merit function equation (3.1) """
>     return dot(x,dot(A,x))-c*dot(x,x)
>
> def meritp(x):
>     """ Gradient of the merit function """
>     return 2*dot(A,x)-2*c*x
>
> def P(y):
>     """ Projection onto \Omega defined as
>     P_\Omega(y) = \argmin\limits_{x \in \Omega} \|x-y\|_2 """
>     n_y = linalg.norm(y)
>     if n_y <=1: return y
>     else: return y/n_y
>
> def dyn(x,t):
>     """ Dynamical system equation (3.2) """
>     tmp = zeros(n,float)
>     tmp = -(x-P(x-meritp(x)))
>     return tmp
>
> n = 6
> L = diag(r_[array(([-2.,-0.2])),zeros(2),ones(n-4)]) # Example 1 p. 46
> B = rand(n,n)
> Q,R = linalg.qr(B)
> A = dot(Q.T,dot(L,Q))
> A = 0.5*(A+A.T)
> c = linalg.norm(A,1)+1. # default value in our numerical computation
> print c
>
> x_0 =-ones(n) # Starting point
> T = linspace(0.,30.,100)
> x = integrate.odeint(dyn,x_0,T)
> plot(T,x[:,0],T,x[:,1],T,x[:,2],T,x[:,3],T,x[:,4],T,x[:,5])
> print dot(x[-1,:],dot(A,x[-1,:])),dot(x[-1,:],x[-1,:])
> show()
>
> ------------------------------------------------------------------------
>
> _______________________________________________
> SciPy-user mailing list
> SciPy-user@scipy.org
> http://projects.scipy.org/mailman/listinfo/scipy-user
>