unitri.py

# Copyright 2020 Jannis Harder
#
# Permission to use, copy, modify, and/or distribute this software for any
# purpose with or without fee is hereby granted.
#
# THE SOFTWARE IS PROVIDED "AS IS" AND THE AUTHOR DISCLAIMS ALL WARRANTIES WITH
# REGARD TO THIS SOFTWARE INCLUDING ALL IMPLIED WARRANTIES OF MERCHANTABILITY
# AND FITNESS. IN NO EVENT SHALL THE AUTHOR BE LIABLE FOR ANY SPECIAL, DIRECT,
# INDIRECT, OR CONSEQUENTIAL DAMAGES OR ANY DAMAGES WHATSOEVER RESULTING FROM
# LOSS OF USE, DATA OR PROFITS, WHETHER IN AN ACTION OF CONTRACT, NEGLIGENCE OR
# OTHER TORTIOUS ACTION, ARISING OUT OF OR IN CONNECTION WITH THE USE OR
# PERFORMANCE OF THIS SOFTWARE.
import numpy as np


def unitri_decomp(M):
    '''
    Unitriangular decomposition.

    Given a non-singular square matrix A find a unitriangular decomposition:

        A = P @ V @ L @ U @ Q

    Parameters
    ----------
    A
        Non-singular square matrix

    Returns
    -------
    P
        Permutation matrix (of dtype int)

    V
        Upper unitriangular matrix, relatively sparse

    L
        Lower unitriangular matrix

    U
        Upper unitriangular matrix, apart from the lowest diagonal element,
        which will be `|det A|` (i.e. unitriangular if `|det A| = 1`).

    Q
        Signed permutation matrix (of dtype int)
    '''
    M = np.copy(M)

    n, m = M.shape
    if n != m:
        raise ValueError('matrix is non-square')

    L = np.identity(n, dtype=M.dtype)
    V = np.identity(n, dtype=M.dtype)

    col_active = np.ones(n, dtype='bool')
    row_active = np.ones(n, dtype='bool')

    row_perm = np.zeros(n, dtype='int')
    col_perm = np.zeros(n, dtype='int')

    for t in range(n - 1):
        W = np.minimum(np.abs(M - 1), np.abs(M + 1))

        W[:, ~col_active] = np.inf
        W[~row_active, :] = np.inf

        i, j = np.unravel_index(np.argmin(W), M.shape)

        target = np.sign(M[i, j])

        row_perm[t] = i
        col_perm[t] = j

        row_active[i] = False
        col_active[j] = False

        if M[i, j] != target:
            k = np.argmax(np.abs(M[:, j]) * row_active)

            a = (target - M[i, j]) / M[k, j]

            M[i, :] += a * M[k, :]
            M[i, j] = target

            V[:, k] -= a * V[:, i]

            L[:, k] -= a * L[:, i]
            L[i, :] += a * L[k, :]

        for l in range(n):
            if row_active[l]:
                b = M[l, j] * target
                M[l, :] -= b * M[i, :]
                M[l, j] = 0

                L[:, i] += b * L[:, l]

    row_perm[n - 1] = np.argmax(row_active)
    col_perm[n - 1] = np.argmax(col_active)

    V = V[row_perm][:, row_perm]
    L = L[row_perm][:, row_perm]
    U = M[row_perm][:, col_perm]

    P = np.identity(n, dtype='int')[:, row_perm]
    Q = np.identity(n, dtype='int')[col_perm]

    for i in range(n):
        if U[i, i] < 0:
            U[:i + 1, i] *= -1
            Q[i, :] *= -1

    return P, V, L, U, Q

Leaving this here so I don't wonder why I have this bookmarked:

There's a graphics trick that used to be widely known that has now probably almost vanished from the graphics consciousness - you can do rotations by applying three shears in a row.

• https://cohost.org/tomforsyth/post/891823-rotation-with-three

it generalizes to higher dimensions and to any volume preserving linear transformation A, not just rotations. So you could use this to linearly transform voxel data or go even higher dimensional. To do that you decompose A (permuting/flipping axes as required) into 3 unitriangular matrices and those are each just a bunch of axis aligned shears (exactly one shear each in the 2d case).

• https://cohost.org/tomforsyth/post/891823-rotation-with-three#comment-0003da45-87c0-411d-a109-aad23ee54ea3