System of bilinear equations look like the following                               y                      T                                    A                      i                          x        =                  g                      i                                   for                     i        =        1        ,        2        ,        …        ,        r                 for some integer                     r                 where                               A                      i                                   are matrices and                               g                      i                                   are some real numbers. These arise in many subjects like engineering, biology, statistics etc.
We consider here the solution theory for bilinear equations in integers. Let the given system of bilinear equation be
                                                                        a                                  x                                      1                                                                    x                                      2                                                  +                b                                  x                                      1                                                                    y                                      2                                                  +                c                                  x                                      2                                                                    y                                      1                                                  +                d                                  y                                      1                                                                    y                                      2                                                                              =                                            α                                                                    e                                  x                                      1                                                                    x                                      2                                                  +                f                                  x                                      1                                                                    y                                      2                                                  +                g                                  x                                      2                                                                    y                                      1                                                  +                h                                  y                                      1                                                                    y                                      2                                                                              =                                            β                                                            This system can be written as
                                          [                                                            a                                                  b                                                  c                                                  d                                                                              e                                                  f                                                  g                                                  h                                                      ]                                                [                                                                                x                                          1                                                                            x                                          2                                                                                                                                        x                                          1                                                                            y                                          2                                                                                                                                        y                                          1                                                                            x                                          2                                                                                                                                        y                                          1                                                                            y                                          2                                                                                            ]                          =                              [                                                            α                                                                              β                                                      ]                                  Once we solve this linear system of equations then by using rank factorization below, we can get a solution for the given bilinear system.
                    m        a        t        (                              [                                                                                x                                          1                                                                            x                                          2                                                                                                                                        x                                          1                                                                            y                                          2                                                                                                                                        y                                          1                                                                            x                                          2                                                                                                                                        y                                          1                                                                            y                                          2                                                                                            ]                          )        =                              [                                                                                x                                          1                                                                            x                                          2                                                                                                            x                                          1                                                                            y                                          2                                                                                                                                        y                                          1                                                                            x                                          2                                                                                                            y                                          1                                                                            y                                          2                                                                                            ]                          =                              [                                                                                x                                          1                                                                                                                                        y                                          1                                                                                            ]                                                [                                                                                x                                          2                                                                                                            y                                          2                                                                                            ]                                  Now we solve first equation by using smith normal form, given any                     m        ×        n                 matrix                     A                , we can get two matrices                     U                 and                     V                 in                                           SL                                m                          (                  Z                )                 and                                           SL                                n                          (                  Z                )                , respectively such that                     U        A        V        =        D                , where                     D                 is as follows:
                    D        =                                            [                                                                                          d                                              1                                                                                                  0                                                        0                                                        …                                                        0                                                                                        0                                                                              d                                              2                                                                                                  0                                                        …                                                        0                                                                                        ⋮                                                                                                                  d                                              s                                                                                                  0                                                                                                          0                                                        0                                                        0                                                        …                                                        0                                                                                        ⋮                                                        ⋮                                                        ⋮                                                        ⋮                                                        ⋮                                                              ]                                            m            ×            n                                  where                               d                      i                          >        0                 and                               d                      i                                    |                          d                      i            +            1                                   for                     i        =        1        ,        2        ,        …        ,        s        −        1                . It is immediate to note that given a system                     A                              x                          =                              b                                  , we can rewrite it as                     D                              y                          =                              c                                  , where                     V                              y                          =                              x                                   and                                           c                          =        U                              b                                  . Solving                     D                              y                          =                              c                                   is easier as the matrix                     D                 is somewhat diagonal. Since we are multiplying with some nonsingular matrices we have the two system of equations to be equivalent in the sense that the solutions of one system have one-to-one correspondence with the solutions of another system. We solve                     D                              y                          =                              c                                  , and take                                           x                          =        V                              y                                  . Let the solution of                     D                              y                          =                              c                                   is
                                          y                          =                              [                                                                                a                                          1                                                                                                                                        b                                          1                                                                                                                    s                                                                              t                                                      ]                                  where                     s        ,        t        ∈                  Z                         are free integers and these are all solutions of                     D                              y                          =                              c                                  . So, any solution of                     A                              x                          =                              b                                   is                     V                              y                                  . Let                     V                 be given by
                    V        =                              [                                                                                a                                          11                                                                                                            a                                          12                                                                                                            a                                          13                                                                                                            a                                          14                                                                                                                                        a                                          21                                                                                                            a                                          22                                                                                                            a                                          23                                                                                                            a                                          24                                                                                                                                        a                                          31                                                                                                            a                                          32                                                                                                            a                                          33                                                                                                            a                                          34                                                                                                                                        a                                          41                                                                                                            a                                          42                                                                                                            a                                          43                                                                                                            a                                          44                                                                                            ]                          =                              [                                                                                A                                          1                                                                                                            B                                          1                                                                                                                                        C                                          1                                                                                                            D                                          1                                                                                            ]                                  Then                                           x                                   is
                    M        =        m        a        t        (                              x                          )        =                              [                                                                                a                                          11                                                                            a                                          1                                                        +                                      a                                          12                                                                            b                                          1                                                        +                                      a                                          13                                                        s                  +                                      a                                          14                                                        t                                                                      a                                          31                                                                            a                                          1                                                        +                                      a                                          32                                                                            b                                          1                                                        +                                      a                                          33                                                        s                  +                                      a                                          34                                                        t                                                                                                  a                                          21                                                                            a                                          1                                                        +                                      a                                          22                                                                            b                                          1                                                        +                                      a                                          23                                                        s                  +                                      a                                          24                                                        t                                                                      a                                          41                                                                            a                                          1                                                        +                                      a                                          42                                                                            b                                          1                                                        +                                      a                                          43                                                        s                  +                                      a                                          44                                                        t                                                      ]                                  We want matrix                     M                 to have rank 1 so that the factorization given in second equation can be done. Solving quadratic equations in 2 variables in integers will give us the solutions for a bilinear systems. This method can be extended to any dimension, but at higher dimension solutions become more complicated. This algorithm can be applied in Sage or Matlab to get to the equations at end.