[SciPy-dev] GSoC weekly report
Alan G Isaac
aisaac@american....
Fri Jul 6 09:10:18 CDT 2007
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
More information about the Scipy-dev
mailing list