bugGNU Octave - Bugs: bug #53140, Solution of a system of linear...


bug #53140: Solution of a system of linear equations takes forever and hurts OS performance.

Submitted by:  None
Submitted on:  Tue 13 Feb 2018 10:48:10 AM UTC  
Category: PerformanceSeverity: 3 - Normal
Priority: 4Item Group: Performance
Status: Ready For TestAssigned to: None
Originator Name: Octave userOriginator Email: -unavailable-
Open/Closed: OpenRelease: 4.2.1
Operating System: Any

Add a New Comment(Rich Markup)

You are not logged in

Please log in, so followups can be emailed to you.


(Jump to the original submission Jump to the original submission)

Wed 21 Feb 2018 11:06:45 AM UTC, comment #17:

@Rik: I set xz to b_nr * b_nc and see

that is a factor 4. The same test in matlab R2017b gives

Therefore, it is no more a problem with memory allocation. And there is a difference between direct-solve and first-factorize, in Matlab too. What I am not sure of, is whether the factorize function in dSparse.cc performs a long [L,U,P,Q,R] factorization or not.

@Rik: you should apply your patch to CSparse.cc, too.

Marco Caliari <caliari>
Project Member
Tue 20 Feb 2018 02:49:33 PM UTC, comment #16:

The original issue, why it was taking 10 hours, was definitely caused by problems of memory allocation.

However, I doubt that this explains the ~2.5X speed advantage of the LU approach. When I had the backslash code instrumented with std::cerr statements I found that the output was resized ~10 times. In each case it didn't seem to take very long.

One way to check would be to hardcode the matrix size to be the final size of x at the beginning of the algorithm. The resize code would never be activated and you could compare just algorithm times.

Rik <rik5>
Project Administrator
Tue 20 Feb 2018 09:26:31 AM UTC, comment #15:

@Rik: is it possible to understand whether it is a problem of memory allocation or not? Because I find quite strange this idea

If you try something like

you will see that x is "never" sparse. I would like to investigate more, but I have no time till the end of February.

Marco Caliari <caliari>
Project Member
Mon 19 Feb 2018 08:00:47 PM UTC, comment #14:

There is probably still an issue about choosing the best algorithm. The method implemented by '\' still seems to be slower than using LU.

Rik <rik5>
Project Administrator
Mon 19 Feb 2018 07:39:15 PM UTC, comment #13:

The problem seems to have been an assumption that floating point arithmetic would be used for mathematical expressions. But because all operands were integers there was no need for the compiler to promote any integer to double. When performing the calculation with integers the result was a value which exceeded the storage size of octave_index_type (signed int), and that resulted in a negative number for the size.

Here is how I recoded it.

I tried to show more clearly that a fraction of the current number of non-zero elements x_nz is added to x_nz to determine the new reserved storage space.

I only cast one of the operands to a double and then relied on the implicit compiler rules. However, if it seems clearer we could also code this using static_casts on every value which we actually want to be double. Or we could use the double constructor explicitly.

I also changed the default increase in size from 10 elements to 100 elements. 10 double values is ~80 bytes which is pretty small given most modern machine's installed memory. Even 100 elements is ~800B or about 1kB which is also pretty small. As an example, I found that the smallest increase in size using the problem matrix was ~8,000 elements.

I checked this change in on the stable branch here (http://hg.savannah.gnu.org/hgweb/octave/rev/956e28867c80).

Marking as Ready for Test.

Rik <rik5>
Project Administrator
Mon 19 Feb 2018 06:19:26 PM UTC, comment #12:

After some more debugging, there is definitely a problem here.

I instrumented the code to be

And when running, I see that the size argument goes negative because octave_idx_type is a sized value, but I can't see why we would ever want to change the capacity of the sparse matrix to a negative value.

Runtime results:

As you can see, the overall size keeps going up by just 10 values which means it takes many iterations before the sparse matrix is increased sufficiently in size.

Rik <rik5>
Project Administrator
Mon 19 Feb 2018 04:41:22 PM UTC, comment #11:

I think it is explainable. jwe has done a lot of work on the symbol table on the development branch. It is possible that calls like change_capacity() do the equivalent of a malloc, rather than a realloc, and therefore Octave starts to run out of memory.

I'll do a little more debugging to see what is actually taking the time.

Rik <rik5>
Project Administrator
Mon 19 Feb 2018 08:35:36 AM UTC, comment #10:


1) this part of code is the same in stable and default branches, and thus comment #2 remains unexplained.

2) in the same file there are other matrix resizing parts, with a slightly different computation of sz. This maybe partially explains comment #4.

Marco Caliari <caliari>
Project Member
Sun 18 Feb 2018 05:51:00 AM UTC, comment #9:

I put in some debug statements and the slowdown is in this code which resizes the matrix

Rik <rik5>
Project Administrator
Fri 16 Feb 2018 07:31:36 AM UTC, comment #8:

@Rik: to_suitesparse_intptr does not exist in stable. What happens on previous versions of octave?

Marco Caliari <caliari>
Project Member
Fri 16 Feb 2018 12:57:26 AM UTC, comment #7:

Not much difference. What happens if you backport just this one function? Does stable then work?

Rik <rik5>
Project Administrator
Thu 15 Feb 2018 08:53:16 AM UTC, comment #6:

If my previous consideration is valid, this is the difference between SparseMatrix::fsolve in stable and default branch

Marco Caliari <caliari>
Project Member
Wed 14 Feb 2018 03:16:12 PM UTC, comment #5:

The problem seems to be the sparsity of b. If I try A \ full (b), then it works in 4.2.1. So, there is something between Matrix SparseMatrix::fsolve and SparseMatrix SparseMatrix::fsolve.

Marco Caliari <caliari>
Project Member
Tue 13 Feb 2018 10:11:10 PM UTC, comment #4:

I changed the test script to be

On my system, the two results are nearly identical (magnitude of max difference is around 1e-20). X has 3 more nonzero entries than XLU, but their absolute values are all around 1e-22.

Performance is significantly different. About 75 seconds for the leftdiv operator and 20 seconds combined for the lu factorization and solve steps. Again, if someone is interested in digging into this, the place to start is SparseMatrix::fsolve.

John W. Eaton <jwe>
Project Administrator
Tue 13 Feb 2018 09:09:33 PM UTC, comment #3:

If someone would like to understand this problem, the place to start is probably the xleftdiv function in libinterp/corefcn/sparse-xdiv.cc. From there, verify that the arguments to the SparseMatrix::solve function are the same in both 4.2.1 and dev. From there, you should end up in the SparseMatrix::fsolve function that takes a SparseMatrix B array. As far as I can tell, there are no significant changes in that function between 4.2.1 and dev.

If the inputs to SparseMatrix::fsolve are the same for both versions, then why does it in 4.2.1 take so much longer than the current one in the dev sources? If the inputs are different, then you'll have to find out why.

John W. Eaton <jwe>
Project Administrator
Tue 13 Feb 2018 06:32:30 PM UTC, comment #2:

Confirmed. This seems to be an issue wit the 4.2.1 release. The development branch that will become 4.4.0 works fine. I tested on a Linux system which proves this is unrelated to the operating system.

I suspect that rather than trying to decipher the cause it may be easier to simply release a new version of Octave which is something that the Maintainer's would like to do anyways.

Rik <rik5>
Project Administrator
Tue 13 Feb 2018 04:53:22 PM UTC, comment #1:

I just tested the matrix with the development version of Octave on Linux:

The lu_test script uses your second code. Clearly this is a big matrix, but I don't see anything like 10 hours.

I'll see if I can start a virtual machine running Windows XP and check the result there.

Rik <rik5>
Project Administrator
Tue 13 Feb 2018 10:48:10 AM UTC, original submission:

The execution of the following code takes at least ten hours or more.

clear all;
close all;

load linsolve_test.mat;

x = A \ b;

However the LU decomposition of A takes only a few seconds:

[L, U, P, Q, R] = lu (A);
x = Q * ( U \ ( L \ ( P * ( R \ b ) ) ) );



(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach File(s):

Attached Files
file #43269:  linsolve_test.7z added by None (1MiB - application/octet-stream)


Depends on the following items: None found

Items that depend on this one: None found


Carbon-Copy List
  • -unavailable- added by caliari (Posted a comment)
  • -unavailable- added by jwe (Posted a comment)
  • -unavailable- added by rik5 (Posted a comment)

    Do you think this task is very important?
    If so, you can click here to add your encouragement to it.
    This task has 0 encouragements so far.

    Only project members can vote.


    Please enter the title of George Orwell's famous dystopian book (it's a date):



    Follow 5 latest changes.

    Date Changed By Updated Field Previous Value => Replaced By
    Mon 19 Feb 2018 07:39:15 PM UTCrik5StatusConfirmed=>Ready For Test
    Tue 13 Feb 2018 10:24:00 PM UTCrik5Priority5 - Normal=>4
    Tue 13 Feb 2018 06:32:30 PM UTCrik5StatusNone=>Confirmed
      Operating SystemMicrosoft Windows=>Any
    Tue 13 Feb 2018 10:48:10 AM UTCNoneAttached File-=>Added linsolve_test.7z, #43269

    Back to the top

    Powered by Savane 3.1-cleanup1