Inverse of a given matrix A of order (n,n).

invmat.py

"""
## LICENSETYPE: Demo
## MODELTYPE: LP


Inverse of a given matrix A of order (n,n).

Two methods are considered.
The first one computes the columns of the inverse x(i) by solving n algebraic
linear systems Ax(i)=e(i), where e(i) are the i-th column of the unity matrix,
i=1,...,n.
The second one determine the inverse by solving the linear system AX=I, where
X is the inverse of A.
"""

from __future__ import annotations

import os

import numpy as np
from gamspy import (
    Alias,
    Container,
    Equation,
    Model,
    Parameter,
    Set,
    Sum,
    Variable,
)


def data_records():
    # a records
    a_recs = np.array(
        [
            [1, 2, 3, 4, 4],
            [1, 3, 4, 3, 1],
            [1, 4, 1, 2, 6],
            [2, 4, 1, 1, 1],
            [3, 1, 5, 2, 7],
        ]
    )

    return a_recs


def main():
    m = Container(
        system_directory=os.getenv("SYSTEM_DIRECTORY", None),
    )

    # SET #
    i = Set(m, name="i", records=[f"i{i}" for i in range(1, 6)])

    # ALIASES #
    j = Alias(m, name="j", alias_with=i)
    k = Alias(m, name="k", alias_with=i)

    # PARAMETERS #
    a = Parameter(
        m,
        name="a",
        domain=[i, j],
        records=data_records(),
        description="matrix to be inverted",
    )
    ident = Parameter(
        m, name="ident", domain=[i, j], description="the identity matrix"
    )

    ident[i, i] = 1

    # Method 1. Solving the systems Ax(i)=e(i)
    # -----------------------------------------

    # VARIABLES #
    obj = Variable(
        m,
        name="obj",
        description="variabile associated to a virtual objective",
    )
    col = Variable(
        m,
        name="col",
        domain=j,
        description="the columns of the inverse matrix",
    )

    # EQUATIONS #
    eobj = Equation(
        m,
        name="eobj",
        type="regular",
        description="name of the virtual objective",
    )
    lin = Equation(
        m,
        name="lin",
        type="regular",
        domain=i,
        description="name of the equations of the systems to be solved",
    )

    # PARAMETERS #
    b = Parameter(m, name="b", domain=i, description="Righ-hand Side term")
    ainv = Parameter(
        m, name="ainv", domain=[i, j], description="inverse matrix of A"
    )

    eobj[...] = obj == 0
    lin[i] = Sum(k, a[i, k] * col[k]) == b[i]

    invmat1 = Model(
        m,
        name="invmat1",
        equations=m.getEquations(),
        problem="lp",
        sense="min",
        objective=obj,
    )

    for jj in j.toList():
        b[i] = ident[i, jj]
        invmat1.solve()
        ainv[i, jj] = col.l[i]

    print("a:  \n", a.pivot().round(3))
    print("ainv:  \n", ainv.pivot().round(3))

    # Checking the inverse
    aainv = Parameter(
        m,
        name="aainv",
        domain=[i, k],
        description="matrix A multiplied by ainv",
    )
    aainv[i, k] = Sum(j, a[i, j] * ainv[j, k])
    print("aainv:  \n", aainv.pivot().round(3))

    # Method 2. Solving the system AX=I
    # ----------------------------------

    # VARIABLE #
    ainverse = Variable(
        m, name="ainverse", domain=[i, j], description="the inverse of A"
    )

    # EQUATION #
    lineq = Equation(m, name="lineq", type="regular", domain=[i, j])

    lineq[i, j] = Sum(k, a[i, k] * ainverse[k, j]) == ident[i, j]

    invmat2 = Model(
        m,
        name="invmat2",
        equations=m.getEquations(),
        problem="lp",
        sense="min",
        objective=obj,
    )
    invmat2.solve()

    print("ainverse:  \n", ainverse.pivot().round(3))

    # Checking the inverse
    aainv[i, k] = Sum(j, a[i, j] * ainv[j, k])
    print("aainv:  \n", aainv.pivot().round(3))

    # End of invmat


if __name__ == "__main__":
    main()