commons-dev mailing list archives

Site index · List index
Message view « Date » · « Thread »
Top « Date » · « Thread »
From Konstantin Berlin <kber...@gmail.com>
Subject Re: [math] refactoring least squares
Date Tue, 25 Feb 2014 14:22:46 GMT
Hi,

I am really having problems believing that matrix copying is the major
problem in an optimization algorithm. Copying is O(N^2) operations. Surely,
for any problem where performance would matter, it is completely dwarfed by
the O(N^3) complexity of actually solving the normal equation.

Also, I think testing should be done on an actual large problem where
scaling issuing would show up. The 1000x2 jaccobian would results in a 2x2
normal equation. Surely this is not a good test case.

Konstantin


On Tue, Feb 25, 2014 at 8:35 AM, Evan Ward <evan.ward@nrl.navy.mil> wrote:

>
> On 02/24/2014 05:47 PM, Gilles wrote:
> > On Mon, 24 Feb 2014 14:41:34 -0500, Evan Ward wrote:
> >> On 02/24/2014 01:23 PM, Gilles wrote:
> >>> On Mon, 24 Feb 2014 11:49:26 -0500, Evan Ward wrote:
> >>>> One way to improve performance would be to provide pre-allocated space
> >>>> for the Jacobian and reuse it for each evaluation.
> >>>
> >>> Do you have actual data to back this statement?
> >>
> >> I did some tests with the CircleVectorial problem from the test cases.
> >> The Jacobian is 1000x2, and I ran it 1000 times until hotspot stopped
> >> making it run faster. The first column is the current state of the code.
> >> The second column is with one less matrix allocation each problem
> >> evaluation.
> >>
> >>        trunk    -1 alloc  %change
> >> lu     0.90 s   0.74      -17%
> >> chol   0.90     0.74      -17%
> >> qr     0.87     0.70      -20%
> >>
> >>
> >> I also see similar reductions in runtime using 1e6 observations and 3
> >> observations. We could save 3-4 allocations per evaluation, which
> >> extrapolates to 60%-80% in run time.
> >
> > I would not have expected such a big impact.
> > How did you set up the benchmark? Actually, it would be a sane starting
> > point to create a "performance" unit test (see e.g.
> > "FastMathTestPerformance"
> > in package "o.a.c.m.util").
>
> The test I used: https://gist.github.com/wardev/3f183dfaaf297dc2462c
>
> I also modified CircleVectorial to use the MultivariateJacobianFunction
> to avoid the translation copies.
>
> >
> >>
> >>>> The
> >>>> LeastSquaresProblem interface would then be:
> >>>>
> >>>> void evaluate(RealVector point, RealVector resultResiduals, RealVector
> >>>> resultJacobian);
> >>>>
> >>>> I'm interested in hearing your ideas on other approaches to solve this
> >>>> issue. Or even if this is an issue worth solving.
> >>>
> >>> Not before we can be sure that in-place modification (rather than
> >>> reallocation) always provides a performance benefit.
> >>
> >> I would like to hear other ideas for improving the performance.
> >
> > Design-wise, it is quite ugly to modify input parameters. I think that it
> > could also hurt in the long run, by preventing other ways to improve
> > performance.
>
> I agree, it smells like fortran. ;)
>
> >
> > Why couldn't the reuse vs reallocate be delegated to implementations of
> > the "MultivariateJacobianFunction" interface?
>
> That is the change I made to produce the performance numbers. It works
> well as long as the application has a single thread. Once it is
> multi-threaded we would have to create a new Function for each thread.
> I'm reluctant to do that because I think not copying the whole problem
> for each thread is a nice feature.
>
> Really I think we need 1 matrix per optimization context. The question
> is: How does the Function, which is several method calls below the
> optimize() method, access the current optimization context? If we reuse
> the matrix at the Function level then we are implicitly tying the
> optimization context to the thread context. By passing the matrices as
> parameters we make the context explicit, but we clutter up the API.
> Perhaps there is a better solution we haven't thought of yet.
>
> > Eventualy, doesn't it boil down to creating "RealVector" and "RealMatrix"
> > implementations that modify and return "this" rather than create a new
> > object?
>
> Interesting idea. I think this will remove the allocation+copy when
> applying weights. That just leaves where matrix decompositions create a
> copy of the input matrix.
>
> Best Regards,
> Evan
>
> >
> >
> > Best Regards,
> > Gilles
> >
> >>>
> >>>> Evan
> >>>
> >>
> >> On 02/24/2014 12:09 PM, Luc Maisonobe wrote:
> >>> Hi Evan,
> >>>
> >>> Le 24/02/2014 17:49, Evan Ward a écrit :
> >>>> I've looked into improving performance further, but it seems any
> >>>> further
> >>>> improvements will need big API changes for memory management.
> >>>>
> >>>> Currently using Gauss-Newton with Cholesky (or LU) requires 4 matrix
> >>>> allocations _each_ evaluation. The objective function initially
> >>>> allocates the Jacobian matrix. Then the weights are applied through
> >>>> matrix multiplication, allocating a new matrix. Computing the normal
> >>>> equations allocates a new matrix to hold the result, and finally the
> >>>> decomposition allocates it's own matrix as a copy. With QR there are
3
> >>>> matrix allocations each model function evaluation, since there is no
> >>>> need to compute the normal equations, but the third allocation+copy
is
> >>>> larger. Some empirical sampling data I've collected with the jvisualvm
> >>>> tool indicates that matrix allocation and copying takes 30% to 80% of
> >>>> the execution time, depending on the dimension of the Jacobian.
> >>>>
> >>>> One way to improve performance would be to provide pre-allocated space
> >>>> for the Jacobian and reuse it for each evaluation. The
> >>>> LeastSquaresProblem interface would then be:
> >>>>
> >>>> void evaluate(RealVector point, RealVector resultResiduals, RealVector
> >>>> resultJacobian);
> >>>>
> >>>> I'm interested in hearing your ideas on other approaches to solve this
> >>>> issue. Or even if this is an issue worth solving.
> >>> Yes, I think this issue is worth solving, especially since we are going
> >>> to ship 3.3 and need to fix as much as possible before the release,
> >>> thus
> >>> avoiding future problems. Everything spotted now is worth fixing now.
> >>>
> >>> Your approach seems reasonable, as long as the work arrays are really
> >>> allocated at the start of the optimization and shared only through a
> >>> few
> >>> documented methods like the one you propose. This would mean we can say
> >>> in the javadoc that these area should be used only to fulfill the API
> >>> requirements and not copied elsewhere, as they *will* be modified as
> >>> the
> >>> algorithm run, and are explicitly devoted to avoid reallocation. I
> >>> guess
> >>> this kind of problems is more important when lots of observations are
> >>> performed, which correspond to very frequent use case (at least in the
> >>> fields I know about).
> >>>
> >>> For the record, what you propose seems similar to what is done in the
> >>> ODE package, as the state vector and its first derivatives are also
> >>> kept
> >>> in preallocated arrays which are reused throughout the integration and
> >>> are used to exchange data between the Apache Commons Math algorithm and
> >>> the user problem to be solved. So it is somehting we already do
> >>> elsewhere.
> >>
> >> OK. We could keep the Evaluation interface, which would just reference
> >> the pre-allocated residuals and matrix. If the result parameters are
> >> null the LSP could allocate a matrix of the correct size automatically.
> >> So then the interface would look like:
> >>
> >> Evaluation evaluate(RealVector point, RealVector resultResiduals,
> >> RealVector resultJacobian);
> >>
> >>>
> >>> best regards,
> >>> Luc
> >>>
> >>>> Best Regards,
> >>>> Evan
> >>>>
> >>>>
> >>>>
> >>>> ---------------------------------------------------------------------
> >>>> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> >>>> For additional commands, e-mail: dev-help@commons.apache.org
> >>>>
> >>>>
> >>>>
> >>>
> >>> ---------------------------------------------------------------------
> >>> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> >>> For additional commands, e-mail: dev-help@commons.apache.org
> >>>
> >>
> >>
> >> On 02/24/2014 01:16 PM, Gilles wrote:
> >>> On Mon, 24 Feb 2014 18:09:26 +0100, Luc Maisonobe wrote:
> >>>> Hi Evan,
> >>>>
> >>>> Le 24/02/2014 17:49, Evan Ward a écrit :
> >>>>> I've looked into improving performance further, but it seems any
> >>>>> further
> >>>>> improvements will need big API changes for memory management.
> >>>>>
> >>>>> Currently using Gauss-Newton with Cholesky (or LU) requires 4 matrix
> >>>>> allocations _each_ evaluation. The objective function initially
> >>>>> allocates the Jacobian matrix. Then the weights are applied through
> >>>>> matrix multiplication, allocating a new matrix. Computing the normal
> >>>>> equations allocates a new matrix to hold the result, and finally
the
> >>>>> decomposition allocates it's own matrix as a copy. With QR there
> >>>>> are 3
> >>>>> matrix allocations each model function evaluation, since there is
no
> >>>>> need to compute the normal equations, but the third
> >>>>> allocation+copy is
> >>>>> larger. Some empirical sampling data I've collected with the
> >>>>> jvisualvm
> >>>>> tool indicates that matrix allocation and copying takes 30% to 80%
of
> >>>>> the execution time, depending on the dimension of the Jacobian.
> >>>>>
> >>>>> One way to improve performance would be to provide pre-allocated
> >>>>> space
> >>>>> for the Jacobian and reuse it for each evaluation. The
> >>>>> LeastSquaresProblem interface would then be:
> >>>>>
> >>>>> void evaluate(RealVector point, RealVector resultResiduals,
> >>>>> RealVector
> >>>>> resultJacobian);
> >>>>>
> >>>>> I'm interested in hearing your ideas on other approaches to solve
> >>>>> this
> >>>>> issue. Or even if this is an issue worth solving.
> >>>>
> >>>> Yes, I think this issue is worth solving, especially since we are
> >>>> going
> >>>> to ship 3.3 and need to fix as much as possible before the release,
> >>>> thus
> >>>> avoiding future problems. Everything spotted now is worth fixing now.
> >>>>
> >>>> Your approach seems reasonable, as long as the work arrays are really
> >>>> allocated at the start of the optimization and shared only through
> >>>> a few
> >>>> documented methods like the one you propose. This would mean we can
> >>>> say
> >>>> in the javadoc that these area should be used only to fulfill the API
> >>>> requirements and not copied elsewhere, as they *will* be modified
> >>>> as the
> >>>> algorithm run, and are explicitly devoted to avoid reallocation. I
> >>>> guess
> >>>> this kind of problems is more important when lots of observations are
> >>>> performed, which correspond to very frequent use case (at least in the
> >>>> fields I know about).
> >>>>
> >>>> For the record, what you propose seems similar to what is done in the
> >>>> ODE package, as the state vector and its first derivatives are also
> >>>> kept
> >>>> in preallocated arrays which are reused throughout the integration and
> >>>> are used to exchange data between the Apache Commons Math algorithm
> >>>> and
> >>>> the user problem to be solved. So it is somehting we already do
> >>>> elsewhere.
> >>>
> >>> If I understand correctly what is being discussed, I do not agree with
> >>> this approach.
> >>>
> >>> The optimization/fitting algorithms must use matrix abstractions.
> >>> If performance improvements can achieved, they must happen at the level
> >>> of the appropriate matrix implementations.
> >>>
> >>
> >> The matrix abstractions will still be used in the interface. As far as I
> >> can tell none of the optimizers or linear algebra classes use the matrix
> >> abstractions internally. For example LU, QR, and Cholesky all copy the
> >> matrix data to an internal double[][]. I tried computing the normal
> >> equation in GaussNewton as j.transpose().multiply(j), but the
> >> performance was bad because j.transpose() creates a copy of the matrix.
> >> That's why we have the current ugly for loop implementation with
> >> getEntry() and setEntry(). Maybe matrix "views" could help solve the
> >> issue.
> >>
> >>>
> >>> Best regards,
> >>> Gilles
> >>>
> >>>
> >>>
> >>> ---------------------------------------------------------------------
> >>> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> >>> For additional commands, e-mail: dev-help@commons.apache.org
> >>>
> >>
> >> Regards,
> >> Evan
> >>
> >>
> >>
> >> ---------------------------------------------------------------------
> >> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> >> For additional commands, e-mail: dev-help@commons.apache.org
> >
> >
> > ---------------------------------------------------------------------
> > To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> > For additional commands, e-mail: dev-help@commons.apache.org
> >
>
>
>
> ---------------------------------------------------------------------
> To unsubscribe, e-mail: dev-unsubscribe@commons.apache.org
> For additional commands, e-mail: dev-help@commons.apache.org
>
>

Mime
  • Unnamed multipart/alternative (inline, None, 0 bytes)
View raw message