#131 - Matrix power via square decomposition, v2 #134
Conversation
| @@ -67,11 +67,12 @@ function IntervalMatrixPower(M::IntervalMatrix{T}) where {T} | |||
| return IntervalMatrixPower(M, M, 1) | |||
| end | |||
|
|
|||
| function IntervalMatrixPower(M::IntervalMatrix{T}, k::Int) where {T} | |||
| function IntervalMatrixPower(M::IntervalMatrix{T}, k::Int; | |||
| algorithm::String="power") where {T} | |||
There was a problem hiding this comment.
Should algorithm be a field of the type?
There was a problem hiding this comment.
I was also thinking about that. It is more general to let you choose at every step, but we could at least have a fallback algorithm defined in the type maybe?
There was a problem hiding this comment.
At first sight it looks better to me if the algorithm is in the type, since it is a property of how to get to M^k, and also how to get higher powers. (Of course i'm assuming that there is no mixed algorithm, but the "mixed" algorithm can be an algorithm by itself).
There was a problem hiding this comment.
Anyway, we can revisit this later.
There was a problem hiding this comment.
(Of course i'm assuming that there is no mixed algorithm, but the "mixed" algorithm can be an algorithm by itself).
Yes and no. If the mixed algorithm follows a strict pattern, it can be added as another algorithm. But adding a new algorithm for any pattern is not sustainable. Furthermore, you can imagine patterns that depend on user knowledge (like "next I need more precision") that you cannot express in a static algorithm. So I conclude that a String as a field is not sufficient. If we make the algorithm a field, it should be its own struct <:IntervalMatrixPowerAlgorithm that you can influence from outside.
| b = k - a^2 | ||
|
|
||
| # recursively compute M^a and M^b | ||
| Mᵏ = square(get(IntervalMatrixPower(pow.M, a; algorithm=algorithm))) |
There was a problem hiding this comment.
Is line 178 the same as get(IntervalMatrixPower(pow.M, a*a) ? (Also, this is the functionality that i was asking in this comment. I'll create an issue).
EDIT: This is now #138.
There was a problem hiding this comment.
My question still is why not use get(IntervalMatrixPower(pow.M, a*a; algorithm=algorithm))) here.
There was a problem hiding this comment.
why not use
get(IntervalMatrixPower(pow.M, a*a; algorithm=algorithm)))here.
Because it is not equivalent. We do not use square for arbitrary square numbers but only for powers of 2.
There was a problem hiding this comment.
Oh, damn, your version was the correct one:
A^(a²) = A^(a*a) = (A^a)^a != (A^a)² for a != 2
Closes #131.
Currently the algorithm for the recursive call cannot be changed.