Notebooks · lib/linalg

An inverse and a reproducible test matrix

Two things that belong under lib/schur, lib/are, lib/lyap, lib/svd and lib/balance rather than inside any one of them.

inverse solves against each column of the identity. It is here for the places where the inverse itself is needed, such as Newton's iteration for the matrix sign function; anywhere the answer is M^-1 b or M^-1 B, call lin_solve on the right-hand side instead and do not form the inverse at all.

test_matrix is a deterministic mess: the same linear congruential generator every time, so a test that passes today passes tomorrow, and entries with no structure for an algorithm to accidentally exploit. Every module in the library tests against matrices from it.

Used with import linalg.

Look inside this program

Atlas
The program’s map: every function and which calls which.
There’s no recorded run of this program yet.

Shift+Enter in a cell runs it and moves to the next. Running a cell runs every cell above it too, as one program.

[1]lines 21–82

#| The inverse of `M`, column by column against the identity.
def inverse(M) {
    return inv(M)
}

#| An `n` by `n` matrix whose entries are a fixed, ugly, reproducible mess in
#| the range -1 to 1. The same `seed` always gives the same matrix.
def test_matrix(n, seed) {
    A = zeros(n, n)
    x = seed
    for i to n {
        for j to n {
            x = mod(1103515245.0 * x + 12345.0, 2147483648.0)
            A[i, j] = x / 1073741824.0 - 1.0
        }
    }
    return A
}

test "the inverse undoes the matrix, both ways round" {
    ncases = 6
    for si to ncases {
        n = 2 + si
        M = test_matrix(n, 7.0 + 13.0 * si)
        for i to n {
            M[i, i] = M[i, i] + 2.5
        }
        Mi = inverse(M)
        assert_true(all(abs(Mi @ M - eye(n)) <= 1.0e-12))
        assert_true(all(abs(M @ Mi - eye(n)) <= 1.0e-12))
    }
}

test "the same seed gives the same matrix, and a different one does not" {
    n = 4
    A = test_matrix(n, 5.0)
    B = test_matrix(n, 5.0)
    C = test_matrix(n, 6.0)
    same = 0.0
    differs = 0.0
    for i to n {
        for j to n {
            same = max(same, abs(A[i, j] - B[i, j]))
            differs = max(differs, abs(A[i, j] - C[i, j]))
        }
    }
    assert_near(same, 0.0, 1.0e-300)
    assert_lt(0.1, differs)
}

test "its entries stay inside the range it promises" {
    n = 9
    A = test_matrix(n, 1234.0)
    for i to n {
        for j to n {
            assert_lt(abs(A[i, j]), 1.0000001)
        }
    }
}

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors

#| The inverse of `M`, column by column against the identity.
def inverse(M) {
    return inv(M)
}

#| An `n` by `n` matrix whose entries are a fixed, ugly, reproducible mess in
#| the range -1 to 1. The same `seed` always gives the same matrix.
def test_matrix(n, seed) {
    A = zeros(n, n)
    x = seed
    for i to n {
        for j to n {
            x = mod(1103515245.0 * x + 12345.0, 2147483648.0)
            A[i, j] = x / 1073741824.0 - 1.0
        }
    }
    return A
}

test "the inverse undoes the matrix, both ways round" {
    ncases = 6
    for si to ncases {
        n = 2 + si
        M = test_matrix(n, 7.0 + 13.0 * si)
        for i to n {
            M[i, i] = M[i, i] + 2.5
        }
        Mi = inverse(M)
        assert_true(all(abs(Mi @ M - eye(n)) <= 1.0e-12))
        assert_true(all(abs(M @ Mi - eye(n)) <= 1.0e-12))
    }
}

test "the same seed gives the same matrix, and a different one does not" {
    n = 4
    A = test_matrix(n, 5.0)
    B = test_matrix(n, 5.0)
    C = test_matrix(n, 6.0)
    same = 0.0
    differs = 0.0
    for i to n {
        for j to n {
            same = max(same, abs(A[i, j] - B[i, j]))
            differs = max(differs, abs(A[i, j] - C[i, j]))
        }
    }
    assert_near(same, 0.0, 1.0e-300)
    assert_lt(0.1, differs)
}

test "its entries stay inside the range it promises" {
    n = 9
    A = test_matrix(n, 1234.0)
    for i to n {
        for j to n {
            assert_lt(abs(A[i, j]), 1.0000001)
        }
    }
}

# SPDX-License-Identifier: MIT OR Apache-2.0
# SPDX-FileCopyrightText: 2026 The Viv authors