[SciPy-dev] GSoC weekly report

[Alan G Isaac]

I hope potential users or contributors will provide Dmitrey 
some feedback on this design decision.

Thank you,
Alan Isaac


On Fri, 06 Jul 2007, dmitrey apparently wrote: 
> 2. Patterns: I should decide which way of implementing those ones is the 
> best. 
> For example, consider the problem 
> n,m,s = 10000,1000,1000 
> ((x+15)^2).sum() -> min# size(x)=n 
> acoording to constraints 
> x1^2 + x2^2 + ... + xm^2 <=c1 
> x2^2 + x3^2 + ... + x[m+1]^2 <=c2 
> ... 
> xs^2+x[s+1]^2+... + x[m+s]^2 <= c[s] 

> so cPattern will look like 
> 1 1 1...1 0 0 0 0 
> 0 1 1 ..1 1 0 0 0 
> 0 0 1 ..1 1 1 0 0 
> ... 
> 0 0 0 ..0 0 0 0 0 
> 0 0 0 ..0 0 0 0 0 

> if I will store dc as 
> 1) dense matrix 
> 2) sparse matrix 
> it will require 
> 1) O(n x s) ~= 10000*1000*sizeof(float) memory bytes 
> 2) O(m x s) = 1000*1000*(sizeof(float) + 2*sizeof(int)) memory bytes 
> However, lots of solvers require just dL/dx = df/dx + dc/dx + dh/dx = 
> grad(f) + sum(dot(mu,c(x))) + sum(dot(lambda,h(x))) 
> so, if dc/dx will be obtained sequentially as dc[i]/dx or dc/dx[i], we 
> can just summarize these ones step by step and amount of memory needed 
> will be O(n) , approximately sizeof(float)*n 
> However, unfortunately not all solvers need just dL/dx, or L is not just 
> f + <mu,c> + <lambda, h>, like for example in Augmented Lagrangian - 
> related solvers. On the other hand, seems like it should always work 
> with fPattern, because f = Fwhole is always just mere sum of f[j] 
> One more problem is how to implement data storing. The most obvious 
> approach is of course double cycle 
> for i in xrange(n) 
>   for j in xrange(len(c))#number of c inequalities 
>     dc[indexes related to the vector (and its' length) returned from 
> c[j](x) evaluation] = (c[j](x)-c0[j])/p.diffInt 
> ('cause c[j] can return several numbers, if we consider the most general 
> case) 
> However, double cycle is too slow for Python. I intend to write single 
> cycle, based on 
> for i in xrange(n) 
>   j = where(cPattern[i]) 
>  C = <somehow evaluate all c[j](x)> 
>  dc[some ind] = (C[some ind] - C0[some ind])/p.diffInt 
> here some ind could be something like [1,2,3] + [4,5] + [8,9,10] + ... + 
> [100, 101] - values obtained from those j: cPattern[j]=1 

> of couse in code I will divide dcmatrix only once, using ufunc /, when 
> whole matrix dc is obtained