[SciPy-User] Unexpectedly large memory usage in scipy.ode class

[Per Nielsen]

Hi all,

I am solving large linear ODE systems using the QuTip python package (
https://code.google.com/p/qutip/) which uses scipy ODE solvers under the
hood. The system is of he form

dydt = L*y,

where L is a large complex sparse matrix, all pretty standard. In this type
of problem the matrix L is the biggest memory user, expected to be much
larger than the solution vector y itself.

Below is the output of a @profile from the memory_profiler package on the
function setting up the ode object, no actual time-stepping is done (the
code can be found here:
https://github.com/qutip/qutip/blob/master/qutip/mesolve.py#L561).

Line #    Mem usage    Increment   Line Contents
================================================
   562                             @profile
   563                             def _mesolve_const(H, rho0, tlist,
c_op_list, expt_ops, args, opt,
   564                                                progress_bar):
   565                                 """!
   566                                 Evolve the density matrix using an
ODE solver, for constant hamiltonian
   567                                 and collapse operators.
   568                                 """
   569    61.961 MB     0.000 MB
   570    61.961 MB     0.000 MB       if debug:
   571                                     print(inspect.stack()[0][3])
   572
   573                                 #
   574                                 # check initial state
   575                                 #
   576    61.961 MB     0.000 MB       if isket(rho0):
   577                                     # if initial state is a ket and
no collapse operator where given,
   578                                     # fallback on the unitary
schrodinger equation solver
   579    61.961 MB     0.000 MB           if len(c_op_list) == 0 and
isoper(H):
   580                                         return _sesolve_const(H,
rho0, tlist, expt_ops, args, opt)
   581
   582                                     # Got a wave function as initial
state: convert to density matrix.
   583    61.973 MB     0.012 MB           rho0 = rho0 * rho0.dag()
   584
   585                                 #
   586                                 # construct liouvillian
   587                                 #
   588    61.973 MB     0.000 MB       if opt.tidy:
   589    61.973 MB     0.000 MB           H = H.tidyup(opt.atol)
   590
   591   327.887 MB   265.914 MB       L = liouvillian_fast(H, c_op_list)
   592
   593                                 #
   594                                 # setup integrator
   595                                 #
   596   343.168 MB    15.281 MB       initial_vector = mat2vec(rho0.full())
   597   343.168 MB     0.000 MB       r = scipy.integrate.ode(cy_ode_rhs)
   598   343.168 MB     0.000 MB       r.set_f_params(L.data.data,
L.data.indices, L.data.indptr)
   599   343.168 MB     0.000 MB       r.set_integrator('zvode',
method=opt.method, order=opt.order,
   600   343.168 MB     0.000 MB                        atol=opt.atol,
rtol=opt.rtol, nsteps=opt.nsteps,
   601   343.168 MB     0.000 MB
 first_step=opt.first_step, min_step=opt.min_step,
   602   343.172 MB     0.004 MB
 max_step=opt.max_step)
   603   572.055 MB   228.883 MB       r.set_initial_value(initial_vector,
tlist[0])
   604
   605                                 #
   606                                 # call generic ODE code
   607                                 #
   608   602.805 MB    30.750 MB       return _generic_ode_solve(r, rho0,
tlist, expt_ops, opt, progress_bar)

On line 591 the L matrix generated and eats a large chunk of memory, as
expected. However, on line 603 setting the initial condition eats an almost
comparable chunk, despite the fact that the initial vector itself only
takes up ~ 15 MB (line 596).

I find this strange, as I would expect that setting the initial condition
would at most increase the memory usage by approximately the size of the
initial vector.

I have tried to reproduce the problem using a minimal script (see
attachment), but here the memory usage is as expected:

Filename: test_ode2.py

Line #    Mem usage    Increment   Line Contents
================================================
     7                             @profile
     8    18.707 MB     0.000 MB   def runode():
     9    18.707 MB     0.000 MB       N = 5000
    10
    11                                 # M = np.random.rand(N, N)
    12   111.230 MB    92.523 MB       M = sparse.rand(N, N, density=0.05,
format='csr') \
    13   198.797 MB    87.566 MB           + 1j * sparse.rand(N, N,
density=0.05, format='csr')
    14   199.031 MB     0.234 MB       y0 = np.random.rand(N, 1) + 1j *
np.random.rand(N, 1)
    15
    16   199.031 MB     0.000 MB       t0 = 0.0
    17
    18   199.031 MB     0.000 MB       def f(t, y, M):
    19                                     # return np.dot(M, y)
    20                                     return M.dot(y)
    21
    22   199.031 MB     0.000 MB       r = ode(f)
    23   199.031 MB     0.000 MB       r.set_integrator('zvode', atol=1e-10)
    24   199.035 MB     0.004 MB       r.set_f_params(M)
    25   199.035 MB     0.000 MB       r.set_initial_value(y0, t0)

Does someone with more insight into the scipy.ode solver might have an idea
of whats going on? I looked in the file myself but didnt not see any
indications of large memory consumptions.

Best,
Per
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20130531/15012e16/attachment.html>
-------------- next part --------------
A non-text attachment was scrubbed...
Name: test_ode2.py
Type: application/octet-stream
Size: 801 bytes
Desc: not available
URL: <http://mail.scipy.org/pipermail/scipy-user/attachments/20130531/15012e16/attachment.obj>