Fast calculation of Nth generalized Fibonacci number of order K?

Yes, there is a faster way to compute this.

Matrix Diagonalisation

The matrix exponentiation can be optimised if the matrix is diagonalisable.

It turns out it is diagonalisable. Indeed, we know that the n-order Fibonacci matrix M is a square matrix and all its rows and columns are linearly independent (by design). This means the rank of the matrix is n (i.e. full-rank matrix). Here is a (more complex) alternative explanation: the Spectral theorem states that a normal matrix can be diagonalised by an orthonormal basis of eigenvectors. Since we can prove that M @ M.T is equal to M.T @ M so M is a normal matrix and so it can be diagonalised.

This means M = Q @ L @ inv(Q) where Q contains the eigenvectors (one per column) and L is a diagonal matrix containing the eigenvalues (one on each item of the diagonal) (see singular value decomposition, a.k.a. SVD, and eigen-decomposition of a matrix for more information). We know that all eigenvalues are non-zeros ones since the matrix is full-rank.

We can compute that in Numpy with L, Q = np.linalg.eig(M). Note that the eigen values and eigenvector can contain complex numbers. We can check that everything so far works as expected with:

D = np.diag(L)
P = np.linalg.inv(Q)
T = Q @ D @ P
np.allclose(M, T)

We can now compute Mⁿ more efficiently than with square exponentiation. Indeed:

Mⁿ = (Q D P)ⁿ
   = (Q D P) (Q D P) (Q D P) ... (Q D P)
   = Q D (P Q) D (P Q) D (P  ...  Q) D P
   = Q D I D I D (I  ...  I) D P                        since P @ Q = I
   = Q D D D ... D P
   = Q Dⁿ P

Q @ Dⁿ @ P can be computed in O(k log n + k³). Indeed, Dⁿ is actually equal to Dn (item-wise exponentiation) in Numpy because D is a diagonal matrix. In fact, we can compute Mⁿ = Q @ Dⁿ @ P with:

Dn = np.diag(Ln)
P = np.linalg.inv(Q)
Mn = Q @ Dn @ P

This method is asymptotically about O(k²) times faster than the O(k³ log n). Actually, for fixed-size floating-point (FP) numbers, it is done in O(k³) time!

Unfortunately, fixed-sized FP numbers are certainly not so accurate for big numbers especially for high-order Fibonacci computation (i.e. performing a SVD computation of large matrices). This method can be seen as a very good approximation for the parameter you use for k < 100. That being said, the computation of Ln fails for values n > 16384 (actually even for much smaller values of n). This is especially true when L contains complex numbers. You can write your own eigen-decomposition of a matrix with arbitrary-large FP numbers (e.g. using the decimal package) though this is tedious to write and the code performance will certainly be far from being optimal (despite the algorithm is efficient). Sympy and mpmath might help for that.

Considering arbitrary-large FP/integer numbers, I expect the computation of the n-order Fibonacci term to be more expensive than expected first. Indeed, AFAIK F(n) is a number having O(n) bits! This means the above method is actually O(n) more expensive with such large numbers assuming basic operations are done in O(n) time. This is the case for the addition but not the multiplication. The native CPython multiplication of large numbers uses the Karatsuba algorithm running in about O(n1.6) time, while decimal.Decimal uses the Number Theoretic Transform algorithm apparently running in O(n log n) time (thanks to @KellyBundy for pointing that out). The optimal complexity for a multiplication of n-bit numbers is O(n log n). As of this sentence, I will assume a O(n log n) algorithm is used for multiplying arbitrary-large numbers (so a code using decimal.Decimal for example). In fact, this also applies to the other algorithm provided in the question: you missed the fact that numbers can be huge!

Considering that, the above algorithm should run in O(k n log² n + n k³) with arbitrary-large FP numbers and considering a O(n log n) multiplication time (for the n-bit numbers).

Improving the diagonalisation-based algorithm

Fortunately, we do not need to compute the whole matrix Mn = Q @ Dn @ P but just a single item. We can replace the last matrix multiplication with a simple dot product. Thus, we only need a column of P and a single line of Q @ Dn. The later can be computed more efficiently since Dn is a diagonal matrix:

Mn[i, j] = np.dot(Q[i]  Ln, P[:,j])
         = np.sum(Q[i]  Ln  P[:,j])
         = np.sum(Q[i]  P[:,j]  Ln)

While we can compute Q[i] P[:,j] efficiently independently of n, we certainly need a very-good precision for this computation since Ln is huge for relatively-big value of n.

This version should still run in O(k log n + k³) for fixed-size numbers (Ln is the expensive part) and it should run in O(k n log² n) for arbitrary-large FP numbers.

I think the complexity of the best algorithm (using arbitrary-large integers) is Ω(k n) since the k-order Fibonacci needs has k terms with a precision of n-bits. This assume the terms are not trivial though. Assuming that, this means this algorithm is actually pretty good already.

Possible further optimisations

One way to improve the diagonalisation-based algorithm further may be to algebraically compute the k-order Fibonacci matrix terms (typically using sympy assuming it can be powerful enough). You might see some terms cancelling each other or some pattern enabling further optimisations. You can find an example here for the classical Fibonacci (with 2x2 matrices). I am not very optimistic since the 2-order algebraic solution found contains 2 irreducible terms, so I expect at least k terms with a k-order Fibonacci. Besides, I think the algorithmic solution found is less numerically stable than the above approach due to a catastrophic cancellation (i.e. ϕⁿ−ψⁿ).