Calculation of consistent initial conditions for ODE

(minor edit)

Hi everyone,

I have written wrapper functions calling MATLAB's ODE solvers. Specifically, my problems are DAEs with large and sparse matrices, similar to matrices coming from finite element analysis. For this reason I use ode15s, ode23t and sometimes ode15i.

When calculating consistent initial conditions for ode15i using decic, everything is fine with small problems, but for even medium problems decic takes ages and may run out of memory. Profiling the code and looking at the implementation of decic I have found that the bottleneck is obviously the QR factorization. However I also have found that the matrix passed to the qr function is converted to full, as shown in the snippet below

function [rank,Q,R,E] = qrank(A)
% QR decomposition with column pivoting and rank determination.
% Have to make A full so that pivoting can be done.
[Q,R,E] = qr(full(A))

The reason of the conversion is explained, however, removing the full function, the code still works.

Is this a legacy option that now is no more required? I don't like to make changes to the built-in files, so in case I would locally copy the dicic function, however, if I remove the full command, am I doing something dangerous?

Below you can find a code that can be used to test decic

% dimensions
nx = 25;
n = 3*nx^2+4*nx+1;

% mass matrix
M = gallery('wathen',nx,nx);

% stiffness matrix
K = gallery('tridiag',n);

% rhs
b = rand(n,1);

% problem: f(t,x,dx) = b(t)-M*dx(t)-K*x(t)
f = @(t,x,dx)b-M*dx-K*x;

% jacobian
options = odeset('Jacobian',{-K,-M});

% consistent initial conditions
x0 = zeros(n,1);
dx0 = zeros(n,1);
[x0,dx0] = decic(f,0,x0,[],dx0,[],options);

% check
norm(f(0,x0,dx0))

[Lateef Adewale Kareem]

Now I understand, the pivoting is needed to ensure that accuracy of the solutin of the linear system is preserved.

What I found is matlab has a function for full matrix and another for sparse matrix. an example is the LDL factoring function. While the full matrx will return all the Ds, and their correcspoing columns of the L, the version that runs on sparse matrix rrequires you to enter the number of Ds you want to be computed.

[Fabio Freschi]

Reading the help of the qr function, the example reported for the 3-output call of the function is done with a sparse matrix. However, it is not explained if the permutaion does something different to pivoting.

[Bjorn Gustavsson]

You can always try. It shouldn't be too complicated to check if the resulting initial condition is valid. When you modify this just remember to rename the function so that you can use either.

[Fabio Freschi]

@Bjorn Gustavsson, thank you for the comment. The check you suggest is already in the code I posted. My point is if the conversion to full matrix is still a necessary step or it is something necessary with old releases of the qr function

[Bjorn Gustavsson]

@Fabio Freschi, yes, but no. You misunderstood what I suggested, what I meant you should check is whether the full-call makes any difference for your initial-condition-calculation. Simply to make your modification of decic without the full in the qrank-function and compare the results with the default decic-function. If there is not and your modification generates OK initial conditions "why bother" (my physicist's general stance to both approximations and rigour). Hopefully I've expressed myself such that it possible to understand this time around?

[Fabio Freschi]

@Bjorn Gustavsson and... yes but no! I actually checked the results with and withoud the full call. My understanding is that the only thing that can be really checked is that the initial conditions are consistent with the problem (my check) because the result is not unique. I also made the check you are suggesting: in all my tests, the results with and without the full call are always identical. However, even if it worked in my limited set of tests, I still bother because I am actually changing a code written by experts without a complete understanding of the implications. And such operations are often a time bomb.

Things would be different if actually the full call refers to old versions of qr, as I believe. decic dates to 2011. Checking the documentation of qr for release 2017a (here), I can find the following

[Q,R,E] = qr(A)

If A is full: [...] The column permutation E is chosen so that abs(diag(R)) is decreasing.

If A is sparse: [...] The column permutation E is chosen to reduce fill-in in R.

There is no such a distinction in the current documentation (here)

[Bjorn Gustavsson]

@Fabio Freschi, OK, now we're at least in the same general region of worrying. If I read decic correctly the use of the QR-output is in lines like these:

  d = - Q' * res;
  dy = E * (R \ d);

Where E is the permutation-matrix. To the best of my understanding the QR-decomposition is a well-defined matrix decomposition, and it would be "very peculiar" if we would get a completely different solution of a linear system of equations if we used

   [Q,R,E] = qr(A);

compared to

    [Q,R,E] = qr(full(A));

The difference in the permutation might/will give normal variations on the level of floating-point accuracy, but surely nothing worse than that (which is very easy for me to say since you're the one working with the time-bomb). If you have a perticularly ill-posed problem that might be a problem, but that would be a difficult problem anyways.

[Fabio Freschi]

Ok, clear. Thank you @Bjorn Gustavsson very much for your time and your help: this is indeed what I wanted to know and be assured of.

The next step should be an official fix from Mathworks, but I don't know how to report it.

I would upvote your comment but it seems it is not possible!

[Bjorn Gustavsson]

Good, this was a learning-step for me too. If you found it useful enough I add my last comment as a reply...

[Christine Tobler]

@Fabio Freschi Quick response on QR and R2017a vs. now:

"Checking the documentation of qr for release 2017a (here), I can find the following

[Q,R,E] = qr(A)

If A is full: [...] The column permutation E is chosen so that abs(diag(R)) is decreasing.

If A is sparse: [...] The column permutation E is chosen to reduce fill-in in R.

There is no such a distinction in the current documentation (here)"

The change in the more recent documentation is purely one of documentation; on rewriting that page in the more modern template, we had preferred making the doc page simpler by removing those additions. I will bring up the possibility of changing this back.

The behavior is still the same as in R2017a: In dense column-permuted QR, you can rely on abs(diag(R)) being decreasing. In the sparse case, that is not the true. However, usually low-rank matrices are detected there in the sense that columns that have norm near some small tolerance are usually dropped completely; this kind of behavior is used in backslash to get an estimated rank of a sparse matrix.

To the larger point of this post, we have received a request to improve performance of decic for the sparse case, and will look into it - thank you for pointing out this issue.

[Fabio Freschi]

@Christine Tobler: Thanks for the comment. It helps in the understanding of the behavior of QR.

The request you mentioned is mine!

Answer 1

[Bjorn Gustavsson]

@Fabio Freschi, OK, now we're at least in the same general region of worrying. If I read decic correctly the use of the QR-output is in lines like these:

  d = - Q' * res;
  dy = E * (R \ d);

Where E is the permutation-matrix. To the best of my understanding the QR-decomposition is a well-defined matrix decomposition, and it would be "very peculiar" if we would get a completely different solution of a linear system of equations if we used

   [Q,R,E] = qr(A);

compared to

    [Q,R,E] = qr(full(A));

The difference in the permutation might/will give normal variations on the level of floating-point accuracy, but surely nothing worse than that (which is very easy for me to say since you're the one working with the time-bomb). If you have a perticularly ill-posed problem that might be a problem, but that would be a difficult problem anyways.