Re: Findminimum too slow for iterative reweighted least
- To: mathgroup at smc.vnet.net
- Subject: [mg123520] Re: Findminimum too slow for iterative reweighted least
- From: DrMajorBob <btreat1 at austin.rr.com>
- Date: Sat, 10 Dec 2011 07:27:26 -0500 (EST)
- Delivered-to: l-mathgroup@mail-archive0.wolfram.com
- References: <201112091055.FAA03827@smc.vnet.net>
- Reply-to: drmajorbob at yahoo.com
The second version, with := not =, takes longer because the Inverse is
computed ALL OVER FROM SCRATCH every time w is needed.
A possible way to speed things up is to simplify, in light of the
following leaf counts:
e.w.e // LeafCount
372180
e.w.e // Simplify // LeafCount
5461
e.w.e // Rationalize // Simplify // LeafCount
891
Count[e.w.e // Rationalize, _Real, Infinity]
0
Similar counts arise for the gradient, so this may improve things:
obj = e.w.e // Rationalize // Simplify;
grad = -2 e.w.j // Rationalize // Simplify;
Timing[irls =
FindMinimum[e.w.e, Sequence @@ startval // Evaluate,
Gradient -> grad]]
{0.486479, {0.165617, {l1 -> 0.550976, l2 -> 1.22172, l3 -> 0.408204,
ps -> 1.42539}}}
compared with
irls = FindMinimum[e.w.e, Sequence @@ startval // Evaluate,
Gradient -> -2 e.w.j]; // Timing
{0.888571, Null}
There's a FindMinimum::lstol error either way, which I haven't thought
much about.
By the way...
I don't like variable names that use the letter l, because it can look TOO
much like the number 1 depending on the font used, how old your eyes are,
etc. The exception is variable names like "list" or "label" or "bill",
where spelling gives us clues.
If you look at the following, quick now... does 11, 12, or 13 appear
anywhere? (The numbers?) It would take a couple of minutes of intense
concentration to decide, so those variable names are a needless source of
confusion.
mat
{{2 (l1^2 + ps)^2, 2 l1 l2 (l1^2 + ps), 2 l1^2 l2^2,
2 l2 l3 (l1^2 + ps), 2 l1 l2^2 l3,
2 l2^2 l3^2}, {2 l1 l2 (l1^2 + ps),
ps (l2^2 + ps) + l1^2 (2 l2^2 + ps), 2 l1 l2 (l2^2 + ps),
l1 l3 (l1^2 + l2^2 + ps), l2 l3 (l1^2 + l2^2 + ps),
2 l1 l2 l3^2}, {2 l1^2 l2^2, 2 l1 l2 (l2^2 + ps), 2 (l2^2 + ps)^2,
2 l1^2 l2 l3, 2 l1 l3 (l2^2 + ps),
2 l1^2 l3^2}, {2 l2 l3 (l1^2 + ps), l1 l3 (l1^2 + l2^2 + ps),
2 l1^2 l2 l3, l2^2 l3^2 + (l1^2 + ps) (l3^2 + ps),
l1 l2 (2 l3^2 + ps), 2 l2 l3 (l3^2 + ps)}, {2 l1 l2^2 l3,
l2 l3 (l1^2 + l2^2 + ps), 2 l1 l3 (l2^2 + ps), l1 l2 (2 l3^2 + ps),
l1^2 l3^2 + (l2^2 + ps) (l3^2 + ps),
2 l1 l3 (l3^2 + ps)}, {2 l2^2 l3^2, 2 l1 l2 l3^2, 2 l1^2 l3^2,
2 l2 l3 (l3^2 + ps), 2 l1 l3 (l3^2 + ps), 2 (l3^2 + ps)^2}}
I can quickly see there are no two digit numbers in this, with or without
resorting to Cases:
mat /. Thread[{l1, l2, l3, ps} -> {m, n, q, p}]
Cases[%, _Integer, Infinity] // Union
{{2 (m^2 + p)^2, 2 m n (m^2 + p), 2 m^2 n^2, 2 n (m^2 + p) q,
2 m n^2 q, 2 n^2 q^2}, {2 m n (m^2 + p),
p (n^2 + p) + m^2 (2 n^2 + p), 2 m n (n^2 + p), m (m^2 + n^2 + p) q,
n (m^2 + n^2 + p) q, 2 m n q^2}, {2 m^2 n^2, 2 m n (n^2 + p),
2 (n^2 + p)^2, 2 m^2 n q, 2 m (n^2 + p) q,
2 m^2 q^2}, {2 n (m^2 + p) q, m (m^2 + n^2 + p) q, 2 m^2 n q,
n^2 q^2 + (m^2 + p) (p + q^2), m n (p + 2 q^2),
2 n q (p + q^2)}, {2 m n^2 q, n (m^2 + n^2 + p) q, 2 m (n^2 + p) q,
m n (p + 2 q^2), m^2 q^2 + (n^2 + p) (p + q^2),
2 m q (p + q^2)}, {2 n^2 q^2, 2 m n q^2, 2 m^2 q^2, 2 n q (p + q^2),
2 m q (p + q^2), 2 (p + q^2)^2}}
{2}
The more VISUALLY simple things are, the more likely it is that you can
come up with ways to further simplify the problem.
Bobby
On Fri, 09 Dec 2011 04:55:51 -0600, Alberto Maydeu-Olivares
<amaydeu at gmail.com> wrote:
> Hi, I'm trying to use FindMinimum to minimize a weighted least squares
> function where the weight matrix needs to be updated at each
> iteration. The weight matrix involves the inverse of a symbolic
> matrix. This should be fast if FindMinimum first evaluates the
> symbolic matrix, then takes the inverse. Given how long it takes, it's
> obviously not doing that. I wonder if anyone can suggest a way to
> force FindMinimum to do that. I obviously tried Evaluate, but does not
> work. At attach a toy example so that you can try. Finding the minimum
> of this toy example takes a full second. It takes minutes to do a
> realistic size job, It should not take more than a couple of seconds.
>
> Thanks a lot for your help.
>
>
> model = {l1^2 + ps, l2*l1, l2^2 + ps, l3*l2, l3*l1, l3^2 + ps};
> theta = {l1, l2, l3, ps};
> j = Outer[D, model, theta];
> data = {2., .42, 3., .56, .48, 1.};
> startval = Transpose[{theta, {.5, .5, .5, .5}}];
> e = data - model;
> mat = {{2 (l1^2 + ps)^2, 2 l1 l2 (l1^2 + ps), 2 l1^2 l2^2,
> 2 l2 l3 (l1^2 + ps), 2 l1 l2^2 l3,
> 2 l2^2 l3^2}, {2 l1 l2 (l1^2 + ps),
> ps (l2^2 + ps) + l1^2 (2 l2^2 + ps), 2 l1 l2 (l2^2 + ps),
> l1 l3 (l1^2 + l2^2 + ps), l2 l3 (l1^2 + l2^2 + ps),
> 2 l1 l2 l3^2}, {2 l1^2 l2^2, 2 l1 l2 (l2^2 + ps), 2 (l2^2 + ps)^2,
> 2 l1^2 l2 l3, 2 l1 l3 (l2^2 + ps),
> 2 l1^2 l3^2}, {2 l2 l3 (l1^2 + ps), l1 l3 (l1^2 + l2^2 + ps),
> 2 l1^2 l2 l3, l2^2 l3^2 + (l1^2 + ps) (l3^2 + ps),
> l1 l2 (2 l3^2 + ps), 2 l2 l3 (l3^2 + ps)}, {2 l1 l2^2 l3,
> l2 l3 (l1^2 + l2^2 + ps), 2 l1 l3 (l2^2 + ps),
> l1 l2 (2 l3^2 + ps), l1^2 l3^2 + (l2^2 + ps) (l3^2 + ps),
> 2 l1 l3 (l3^2 + ps)}, {2 l2^2 l3^2, 2 l1 l2 l3^2, 2 l1^2 l3^2,
> 2 l2 l3 (l3^2 + ps), 2 l1 l3 (l3^2 + ps), 2 (l3^2 + ps)^2}};
>
> (*brute force approach to benchmark*)
> w = Inverse[mat];
> irls = FindMinimum[e. w . e, Evaluate[Sequence @@ startval], Gradient
> -> -2 e. w .j]; // Timing
>
> (*this should work, in fact, it takes almost twice as long*)
> w := Inverse[mat];
> irls = FindMinimum[Evaluate[e. w . e], Evaluate[Sequence @@
> startval], Gradient -> Evaluate[-2 e. w .j]]; // Timing
>
--
DrMajorBob at yahoo.com
- References:
- Findminimum too slow for iterative reweighted least squares
- From: Alberto Maydeu-Olivares <amaydeu@gmail.com>
- Findminimum too slow for iterative reweighted least squares