Handle large integer exponents in matrix powers

[jishnub]

The idea here is that A^p is evaluated as A^Int(p) if isinteger(p), but Int(p) may error if p is large. In such cases, we may express the exponent as p = m*q + r, where m and r are small enough so that Int(m) and Int(r) succeed, and the operation becomes (A^m)^q * A^r. This is evaluated recursively by repeatedly dividing the exponent.

After this, the following works:

julia> [1 1e-10; 0 1] ^ (1e20)
2×2 Matrix{Float64}:
 1.0  1.0e10
 0.0  1.0

Fixes JuliaLang/LinearAlgebra.jl#1082, but I have not carried out any error analysis to check the accuracy of results in general. Suggestions are welcome. In particular, I'm unsure how the rounding in divrem impacts the results if the exponents are huge.

[developfast]

overall lgtm, just suggested modifying some cases

[developfast]

perhaps explicitly adding a case for p==0 to return the identity matrix

[jishnub]

I think this is handled by the A^0 call, which would also ensure type-stability.

[developfast]

hmm ok fair

[developfast]

nit: this one can directly be q>0

[mikmoore]

But the ==1 case is a no-op, so nicer to skip it.

[mikmoore]

Suggested change

	m = typemax(Int)
	m = prevpow(2, typemax(UInt))

typemax(Int) will not be accurate here (although arguably is not that far off). Note the inaccuracy of divrem(2.0^64, typemax(Int)) == (2.0, 0.0) due to the promotion Int -> Float64. Also, I believe the cost of power_by_squaring scales partly with the count_ones in the power.

[longemen3000]

I think that the correct value here is floatintmax(p)

[mikmoore]

I think you're right that it needs to be <=maxintfloat(typeof(p)), otherwise the rem could be wrong. So I think it needs to be UInt(min(prevpow(2, typemax(UInt)), maxintfloat(typeof(p)))).

[developfast]

For the divrem issue, suppose p = 5*10^18 and m = typemax(Int) (which, for a 64-bit system, is 2^63 - 1 or approximately 9.22 × 10^18).

The result of divrem(p, m) would yield q = 0 and r = 5*10^18 (since p < m).
The computation effectively becomes A^(5*10^18), which is not broken down further because r still represents a large exponent. We could choose a smaller m but the tradeoff is using more recursion calls.

[mikmoore]

The flow of all of this seems a little more complicated than necessary (maybe I'm wrong). What about just computing q, r = divrem(p, m) initially and branching on q > 0 rather than p <= m? It seems that could remove at least one if block.

[mikmoore]

Suggested change

	_integerpow(A::AbstractMatrix, p) = A^Integer(p)
	function _integerpow(A::AbstractMatrix, p::Union{Float32, Float64})
	_integerpow(A::AbstractMatrix, p::Integer) = A^p
	function _integerpow(A::AbstractMatrix, p)

It seems that the long definition here is the more general? Apart from Float16 (although even there the overhead should be tiny), I'd feel safer with generic types avoiding the assumed-valid Integer conversion since that's how we got here in the first place.

EDIT: probably my proposed _integerpow(A::AbstractMatrix, p::Integer) is redundant altogether, if not directly overwriting another definition. I haven't traced the whole dispatch logic here.

[stevengj]

Wouldn't a much simpler solution just be to replace isinteger(p) && return integerpow(A, p) with

isinteger(p) && typemin(Int) ≤ p ≤ typemax(Int) && return integerpow(A, p)

? i.e. just treat p the same as any floating-point exponent if it's too big?

[mikmoore]

just treat p the same as any floating-point exponent if it's too big?

What about if the matrix is the negative identity, a permutation matrix, or a similarly "simple" matrix? Using the floating point exponentiation algorithm might cause a poor result. Although I can't say for sure because I haven't looked at the generic algorithm. Is it eigendecomposition based? I guess that would probably work for matrices with trivial diagonalizations (like my examples), but perhaps would cause some undesirable roundoff and error accumulation in others.

[stevengj]

Eigendecompositions are only used for Hermitian matrices; for other matrices it's too dangerous because they might be close to defective.

But nevermind, even schurpow needs to call integerpow, so you can't replace the latter with the former.

Note that this PR needs to move to https://github.com/JuliaLang/LinearAlgebra.jl, I think?