Last active
August 26, 2018 10:47
-
-
Save dpiponi/8983fb46ff84fdb5fc53129ac166aa7a to your computer and use it in GitHub Desktop.
Invert infinite upper triangular matrix
This file contains bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
| > scale :: Num a => a -> [a] -> [a] | |
| > scale a bs = map (a *) bs | |
| Invert an infinite upper triangular matrix with 1 on diagonal | |
| [ | |
| [1, b₀₁, b₀₂, b₀₃, …], | |
| [1, b₁₂, b₁₃, …], | |
| [1, b₂₃, …], | |
| … | |
| ] | |
| giving upper triangular matrix in same format. | |
| > invertUT1 :: Num a => [[a]] -> [[a]] | |
| > invertUT1 ((b00 : b01) : | |
| > b11) = (1 : staggered_sum (zipWith scale (map negate b01) a11)) : | |
| > a11 | |
| The bottom right, `b11` of matrix is inverted recursively in `a11` and this becomes | |
| the bottom right of result. But it's also used to subtract elements from top row | |
| of result. | |
| There is of course no "base case" and yet the recursion works. We also use | |
| some of `a11` to compute the top row. | |
| > where a11 = invertUT1 b11 | |
| > staggered_sum ((c00 : c01) : | |
| > c11) = c00 : zipWith (+) c01 (staggered_sum c11) | |
| Invert an infinite upper triangular matrix | |
| [ | |
| [b₀₀, b₀₁, b₀₂, b₀₃, …], | |
| [b₁₁, b₁₂, b₁₃, …], | |
| [b₂₂, b₂₃, …], | |
| … | |
| ] | |
| Prescales rows by bᵢᵢ⁻¹, applies `invertUT1`, and then postscales columns by bᵢᵢ. | |
| > invertUT :: Fractional a => [[a]] -> [[a]] | |
| > invertUT b = postscale rdiags $ invertUT1 $ zipWith scale rdiags b | |
| > where rdiags = map (recip . head) b | |
| > postscale s@(_ : ss) (b : bs) = zipWith (*) s b : postscale ss bs | |
| "Signed" Stirling numbers of first kind | |
| > stirling1 0 0 = 1 | |
| > stirling1 _ 0 = 0 | |
| > stirling1 0 _ = 0 | |
| > stirling1 n k = -(n-1)*stirling1 (n-1) k+stirling1 (n-1) (k-1) | |
| Stirling numbers of second kind | |
| > stirling2 0 0 = 1 | |
| > stirling2 _ 0 = 0 | |
| > stirling2 0 _ = 0 | |
| > stirling2 n k = k*stirling2 (n-1) k+stirling2 (n-1) (k-1) | |
| > m = [[stirling1 n k | n <- [k..]] | k <- [0..]] | |
| > m' = [[stirling2 n k | n <- [k..]] | k <- [0..]] | |
| > top_left :: Int -> [[a]] -> [[a]] | |
| > top_left n b = map (take n) (take n b) | |
| The two matrices of Stirling numbers are mutual inverses | |
| > main = do | |
| > print $ top_left 6 m | |
| > print $ top_left 6 m' | |
| > print $ top_left 6 $ invertUT1 m |
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment