gauss

det_gauss(A, progress=set())

Return the determinant.

\[ \det A \qquad \mathbb{K}^{N\times N}\to\mathbb{K} \]

Uses Gaussian elimination with complete pivoting.

The matrix will be transformed in-place into an upper triangular matrix (columns left of pivot won't be reduced).

Complexity

There will be at most

• \(\frac{N(N^2-1)}{3}\) scalar subtractions (-),
• \(\begin{cases} \frac{N(N^2+2)}{3}-1 & N>0 \\ 0 & N=0 \end{cases}\) scalar multiplications (*) &
• \(\frac{N(N-1)}{2}\) scalar true divisions (/).

References

Wikipedia - Gaussian elimination - Computing determinants

Source code in linalg\gauss.py

def det_gauss(A:np.ndarray, progress:set[str]|Progress=set()) -> Any:
    r"""Return the determinant.

    $$
        \det A \qquad \mathbb{K}^{N\times N}\to\mathbb{K}
    $$

    Uses Gaussian elimination with complete pivoting.

    The matrix will be transformed in-place into an upper triangular matrix
    (columns left of pivot won't be reduced).

    Complexity
    ----------
    There will be at most

    - $\frac{N(N^2-1)}{3}$ scalar subtractions (`-`),
    - $\begin{cases} \frac{N(N^2+2)}{3}-1 & N>0 \\ 0 & N=0 \end{cases}$ scalar multiplications (`*`) &
    - $\frac{N(N-1)}{2}$ scalar true divisions (`/`).

    References
    ----------
    [Wikipedia - Gaussian elimination - Computing determinants](https://en.wikipedia.org/wiki/Gaussian_elimination#Computing_determinants)
    """
    assert_sqmatrix(A)

    def _det_gauss(A:np.ndarray, progress:Progress) -> Any:
        s:bool = True
        for i in range(A.shape[0]):
            #pivot
            i_max, j_max = \
                    np.unravel_index(np.argmax(np.abs(A[i:, i:])), A[i:, i:].shape)
            if not A[i+i_max, i+j_max]:
                return +A[i+i_max, i+j_max]
            swap_pivot(A, i, i+i_max, i+j_max)
            s ^= bool(i_max) != bool(j_max)
            #reduce (not left of pivot, these elements will not influence result)
            #A[i+1:, i:] -= A[i, i:] * (A[i+1:, i] / A[i, i])[:, np.newaxis]
            vecisub(A[i+1:, i:],
                outer(vectruediv(A[i+1:, i], A[i, i], progress),
                    A[i, i:], progress), progress)
        return posneg(progress.prod_default(np.diag(A)), s)

    if isinstance(progress, set):
        N:int = A.shape[0]
        totals:dict[str,int] = {
            '-': N*(N**2-1)//3,
            '*': N*(N**2+2)//3-1,
            '/': N*(N-1)//2
        }
        with Progress({o:totals[o] for o in progress}, 'det_gauss ') as progress:
            return _det_gauss(A, progress)
    elif isinstance(progress, Progress):
        return _det_gauss(A, progress)
    else:
        raise TypeError('expected progress as `set` or `Progress`')

inv_gauss(A, progress=set())

Return the inverse.

\[ A^{-1} \qquad \mathbb{K}^{N\times N}\to\mathbb{K}^{N\times N} \]

Uses Gaussian elimination with complete pivoting. The matrix will be transformed in-place into the identity matrix.

TODO: Still unnecessary operations like A[i, i]/A[i, i].

Complexity

There here will be

• \(2N^3-2N^2\) scalar subtractions (-),
• \(2N^3-2N^2\) scalar multiplications (*),
• \(2N^2\) scalar true divisions (/),

so \(4N^3-2N^2\) arithmetic operations in total.

References

StackExchange - Gauss-Jordan: Effect of column pivoting on result matrix

Source code in linalg\gauss.py

def inv_gauss(A:np.ndarray, progress:set|Progress=set()) -> np.ndarray:
    r"""Return the inverse.

    $$
        A^{-1} \qquad \mathbb{K}^{N\times N}\to\mathbb{K}^{N\times N}
    $$

    Uses Gaussian elimination with complete pivoting.
    The matrix will be transformed in-place into the identity matrix.

    TODO: Still unnecessary operations like `A[i, i]/A[i, i]`.

    Complexity
    ----------
    There here will be

    - $2N^3-2N^2$ scalar subtractions (`-`),
    - $2N^3-2N^2$ scalar multiplications (`*`),
    - $2N^2$ scalar true divisions (`/`),

    so $4N^3-2N^2$ arithmetic operations in total.

    References
    ----------
    [StackExchange - Gauss-Jordan: Effect of column pivoting on result matrix](https://math.stackexchange.com/a/744213/1170417)
    """
    assert_sqmatrix(A)

    def _inv_gauss(A:np.ndarray, progress:Progress) -> np.ndarray:
        P:np.ndarray = np.eye(A.shape[0], dtype=A.dtype)
        Q:list[int] = list(range(A.shape[0]))

        for i in range(A.shape[0]):
            #pivot
            i_max, j_max = \
                    np.unravel_index(np.argmax(np.abs(A[i:, i:])), A[i:, i:].shape)
            if not A[i+i_max, i+j_max]:
                raise ZeroDivisionError
            swap_pivot(A, i, i+i_max, i+j_max)
            swap_rows(P, i, i+i_max)#, swap_columns(Q, i, i+j_max)
            Q[i], Q[i+j_max] = Q[i+j_max], Q[i]

            #normalize pivot
            #P[i, :] /= A[i, i]
            #A[i, :] /= A[i, i]
            vecitruediv(P[i, :], A[i, i], progress)
            vecitruediv(A[i, :], A[i, i], progress)

            #zeros above and below
            #P[:i, :] -= P[i, :] * A[:i, i][:, np.newaxis]
            #A[:i, :] -= A[i, :] * A[:i, i][:, np.newaxis]
            #P[i+1:, :] -= P[i, :] * A[i+1:, i][:, np.newaxis]
            #A[i+1:, :] -= A[i, :] * A[i+1:, i][:, np.newaxis]
            vecisub(P[:i, :], outer(A[:i, i], P[i, :], progress), progress)
            vecisub(A[:i, :], outer(A[:i, i], A[i, :], progress), progress)
            vecisub(P[i+1:, :], outer(A[i+1:, i], P[i, :], progress), progress)
            vecisub(A[i+1:, :], outer(A[i+1:, i], A[i, :], progress), progress)

        #return matmul(Q, P, progress)
        return P[np.argsort(Q),:]

    if isinstance(progress, set):
        N:int = A.shape[0]
        totals:dict[str,int] = {
            '-': N**2*(2*N-2),
            '*': N**2*(2*N-2),
            '/': 2*N**2
        }
        with Progress({o:totals[o] for o in progress}, 'det_gauss ') as progress:
            return _inv_gauss(A, progress)
    elif isinstance(progress, Progress):
        return _inv_gauss(A, progress)
    else:
        raise TypeError('expected progress as `set` or `Progress`')