BinaryPowerMethod should cache M^**(2^k) steps

[mahrud]

Using BinaryPowerMethod for tensor powers is pretty good, but if one needs every tensor power from 1 to d for instance, then there's a lot of repeated computation of M^**(2^k) for k < log_2 d. I think some opportunistic caching here might go a long way. Not sure exactly how to do it at the generality of BinaryPowerMethod however.

[pzinn]

In LieTypes I actually disabled BinaryPowerMethod:

-- the implementation below seems more reasonable but it's actually slower in most circumstances
LieAlgebraModule ^** ZZ := LieAlgebraModule => (M, n) -> BinaryPowerMethod(M, n, tensor,
    M -> trivialModule M#"LieAlgebra",
    M -> error "LieAlgebraModule ^** ZZ: expected non-negative integer")

after finding out that it was slower than the stupid method:

LieAlgebraModule ^** ZZ := (M, n) -> M.cache#(symbol ^**, n) ??= (
        if n<0 then "error nonnegative powers only";
        if n==0 then trivialModule M#"LieAlgebra"
        else if n==1 then M
        else M**(M^**(n-1)) -- order matters for speed purposes
    )

I'm not sure I understand why exactly but I'm worried there might be other circumstances where the same is true.

[pzinn]

Incidentally, I also don't understand why the last 2 arguments of BinaryPowerMethod are functions.

[mahrud]

Can you give an example where it's slower? My guess it that your M ** (M ^** (n-1)) implementation is faster on the second attempt because of caching.

Incidentally, I also don't understand why the last 2 arguments of BinaryPowerMethod are functions.

As opposed to what? They need some information about the input, so they have to be functions.

[pzinn]

needsPackage"LieTypes"
g=𝔣_4
L=adjointModule g;
time(L^**7;)
-- dumb method:       used 5.24386s (cpu); 3.20123s (thread); 0s (gc)
-- BinaryPowerMethod: used 47.3824s (cpu); 26.455s (thread); 0s (gc)
restart
time(L^**8;)
-- dumb method:       used 7.21557s (cpu); 4.77082s (thread); 0s (gc)
-- BinaryPowerMethod: used 447.432s (cpu); 129.426s (thread); 0s (gc)
needsPackage"LieTypes"
g=𝔞_2
V=LL_(5,3)g
time(V^**15;)
-- dumb method:       used 27.238s (cpu); 17.5987s (thread); 0s (gc)
-- BinaryPowerMethod: used 39.4362s (cpu); 25.7623s (thread); 0s (gc)

[mahrud]

In your case this is because the iterative algorithm only computes characters of L because of the way tensor is implemented, whereas the binary algorithm computes characters for L^**2 and L^**4 and so on as well, which is the heaviest step. However, if I'm not mistaken the character is multiplicative, so a simple precomputation should close most of the gap:

@@ -1205,7 +1207,11 @@ LieAlgebraModule ** LieAlgebraModule := (V,W) -> ( -- cf Humpheys' intro to LA &
 	    	    );
 		if i === null then add(u-rho,a*b*t);
 		)));
-    new LieAlgebraModule from (g,ans)
+    M := new LieAlgebraModule from (g,ans);
+    if V.cache#?character and W.cache#?character
+    then M.cache#character = character V * character W;
+    if V === W then V.cache#(symbol ^**, 2) = M;
+    M

After this, the next bottleneck is this step:

M2/M2/Macaulay2/packages/LieTypes.m2

Lines 1198 to 1199 in 00b89ec

	scanPairs(W#"DecompositionIntoIrreducibles", (w,a) -> -- loop over highest weights of W
	scanPairs(wd, (v,b) -> ( -- loop over all weights of V

In the iterative algorithm, since wd is the weight diagram of the smaller module, there are fewer weights here to deal with, but in the binary algorithm there are potentially many weights to deal with here. I'm not sure what's the algorithm doing here. Is it not possible to loop only over the highest weights for both V and W?

Beyond this it's a question of what you want to optimize:

• the iterative algorithm hogs memory by caching everything from L^**2 .. L^**n without pre-computing characters, but then L^**(n+1) is very fast;
• the binary algorithm only needs to cache L^**(2^k) but then has to compute log(n) many not so easy tensors to give the final result.

Given how many times you use memoize, it seems you don't care about memory, so probably the binary method doesn't make sense for you. However, for general modules the same issues don't apply.

[pzinn]

Is it not possible to loop only over the highest weights for both V and W?

no that'd be too easy...

On the other hand, we typically don't want to compute the characters of V and W when doing V ** W, because these characters can get very big. In fact that would give a trivial algorithm to compute V ** W (subtract inductively from ch(V) ch(W) irreducible characters), but that's not the algorithm we're using.