# On the classification of skew Hadamard matrices of order $\boldsymbol{36}$
and related structures

Makoto Araya${}^{{}^{1}}$ Masaaki Harada${}^{{}^{2}}$ Hadi
Kharaghani${}^{{}^{3}}$
Ali Mohammadian${}^{{}^{4}}$ Behruz Tayfeh-Rezaie${}^{{}^{5}}$
${}^{{}^{1}}$Department of Computer Science, Shizuoka University,
Hamamatsu 432-8011, Japan
${}^{{}^{2}}$Research Center for Pure and Applied Mathematics,
Graduate School of Information Sciences, Tohoku University,
Sendai 980-8579, Japan
${}^{{}^{3}}$Department of Mathematics and Computer Science, University of
Lethbridge,
Lethbridge, Alberta, T1K 3M4, Canada
${}^{{}^{4}}$School of Mathematical Sciences, Anhui University,
Hefei 230601, Anhui, China
${}^{{}^{5}}$School of Mathematics, Institute for Research in Fundamental
Sciences (IPM),
P.O. Box 19395-5746, Tehran, Iran
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
<EMAIL_ADDRESS><EMAIL_ADDRESS>
###### Abstract

Two skew Hadamard matrices are considered SH-equivalent if they are similar by
a signed permutation matrix. This paper determines the number of SH-
inequivalent skew Hadamard matrices of order $36$ for some types. We also
study ternary self-dual codes and association schemes constructed from the
skew Hadamard matrices of order $36$.
Keywords: Association Scheme, Classification, Skew Hadamard Matrix, Ternary
Self-Dual Code.
AMS 2020 Mathematics Subject Classification: 05B20, 15B34, 94B05.

##  1 Introduction

A Hadamard matrix of order $n$ is an $n\times n$ matrix $H$ with entries in
$\\{-1,1\\}$ such that $HH^{\top}=nI$, where $H^{\top}$ is the transpose of
$H$ and $I$ is the identity matrix. It is well known that the order of a
Hadamard matrix is $1$, $2$, or a multiple of $4$ [20]. It is a longstanding
folklore conjecture that Hadamard matrices of order $n$ exist for any $n$
divisible by $4$. The smallest order for which no known Hadamard matrix is
$668$ [13].

A Hadamard matrix $H$ is said to be skew Hadamard if $H+H^{\top}=2I$. Skew
Hadamard matrices are equivalent to doubly regular tournaments [3]. They are
used to construct several combinatorial objects, such as association schemes,
self-dual codes, strongly regular graphs, and more. It is conjectured that
skew Hadamard matrices of order $n$ exist for any $n$ divisible by $4$ [23].
The smallest unknown order of a skew Hadamard matrix is $356$ [5].

A permutation matrix (respectively, signed permutation matrix) is a matrix
with entries in $\\{0,1\\}$ (respectively, $\\{-1,0,1\\}$) which has exactly
one nonzero entry in each row and each column. Two matrices $X$ and $Y$ with
entries in $\\{-1,1\\}$ are said to be H-equivalent if there exist two signed
permutation matrices $P$ and $Q$ such that $Y=PXQ$, otherwise, they are called
H-inequivalent. Two matrices $X$ and $Y$ with entries in $\\{-1,1\\}$ are said
to be SH-equivalent if there exists a signed permutation matrix $P$ such that
$Y=P^{\top}XP$, otherwise, they are called SH-inequivalent.

All H-inequivalent Hadamard matrices of orders up to $32$ have been classified
[11]. Also, all H-inequivalent and SH-inequivalent skew Hadamard matrices of
orders up to $32$ have been classified [8]. The resulting classification is
shown in Table 1. The complete classification of Hadamard matrices of order
$36$ seems to be difficult and perhaps inaccessible.

order | 1 | 2 | 4 | 8 | 12 | 16 | 20 | 24 | 28 | 32
---|---|---|---|---|---|---|---|---|---|---
# H-inequivalent Hadamard matrices | 1 | 1 | 1 | 1 | 1 | 5 | 3 | 60 | 487 | 13710027
# H-inequivalent skew Hadamard matrices | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 16 | 54 | 6662
# SH-inequivalent skew Hadamard matrices | 1 | 1 | 1 | 1 | 1 | 2 | 2 | 16 | 65 | 7227
Table 1: The number of H-inequivalent Hadamard matrices, H-inequivalent skew
Hadamard matrices, and SH-inequivalent skew Hadamard matrices of orders up to
$32$.

This paper determines the number of SH-inequivalent skew Hadamard matrices of
order $36$ for some types. As an application, we study ternary self-dual codes
and association schemes constructed from the skew Hadamard matrices of order
$36$.

##  2 Preliminaries

In this section, we fix our notation and present some preliminary results.
Throughout the paper, we denote all one vector of length $r$ by
$\mathbb{1}_{r}$. Let $H=[h_{uv}]$ be a Hadamard matrix of order $n$. We know
from [4] that, by a sequence of column negations and column permutations,
every four distinct rows $i,j,k,\ell$ of $H$ may be transformed to the form

$\displaystyle\begin{array}[]{rrrrrrrrrr}\mathbbm{1}_{s}&\mathbbm{1}_{t}&\mathbbm{1}_{t}&\mathbbm{1}_{s}&\mathbbm{1}_{t}&\mathbbm{1}_{s}&\mathbbm{1}_{s}&\mathbbm{1}_{t}\\\
\mathbbm{1}_{s}&\mathbbm{1}_{t}&\mathbbm{1}_{t}&\mathbbm{1}_{s}&-\mathbbm{1}_{t}&-\mathbbm{1}_{s}&-\mathbbm{1}_{s}&-\mathbbm{1}_{t}\\\
\mathbbm{1}_{s}&\mathbbm{1}_{t}&-\mathbbm{1}_{t}&-\mathbbm{1}_{s}&\mathbbm{1}_{t}&\mathbbm{1}_{s}&-\mathbbm{1}_{s}&-\mathbbm{1}_{t}\\\
\mathbbm{1}_{s}&-\mathbbm{1}_{t}&\mathbbm{1}_{t}&-\mathbbm{1}_{s}&\mathbbm{1}_{t}&-\mathbbm{1}_{s}&\mathbbm{1}_{s}&-\mathbbm{1}_{t}\end{array}$
(5)

for some uniquely determined $s,t$ with $s+t=\frac{n}{4}$. By negation of the
last row in (5) and then a suitable column permutation if necessary, we may
assume that $s\geqslant t$ and so $0\leqslant
t\leqslant\lfloor\frac{n}{8}\rfloor$. Following [14], we define the type of
the four rows $i,j,k,\ell$ as $T_{ijk\ell}=t$. It is straightforward to check
that $T_{ijk\ell}=\frac{n-|P_{ijk\ell}|}{8}$, where

$P_{ijk\ell}=\sum_{r=1}^{n}h_{ir}h_{jr}h_{kr}h_{\ell r}.$

This shows that ‘type’ is an equivalence invariant, meaning that any negation
of rows and columns and any permutation of columns leaves the type unchanged.
As $H^{\top}$ is a Hadamard matrix, we may define ‘type’ for any quadruple of
columns of $H$ in a similar way. A Hadamard matrix is of type $t$ if it has a
quadruple of rows of type $t$ and no quadruple of rows of type less than $t$.
We refer to [19] for more information about types of Hadamard matrices.

Below, we find a nice form for a quadruple of rows of a skew Hadamard matrix.
It is not hard to verify that any $4\times 4$ skew-symmetric matrix with
entries in $\\{-1,1\\}$ and with constant diagonal is H-equivalent to one of
the matrices

$A=\left[\begin{array}[]{rrrr}1&1&1&1\\\ -1&1&1&1\\\ -1&-1&1&1\\\
-1&-1&-1&1\end{array}\right]\quad\text{ or }\quad
B=\left[\begin{array}[]{rrrr}1&1&1&1\\\ -1&1&1&-1\\\ -1&-1&1&1\\\
-1&1&-1&1\end{array}\right].$ (6)

Therefore, by a sequence of row negations and row permutations, and
simultaneously, column negations and column permutations, every quadruple of
rows of a skew Hadamard matrix may be transformed to one of the forms

$\displaystyle\begin{array}[]{rrrrrrrrrrrrrrrrrr}1&1&1&1&\mathbbm{1}_{s-1}&\mathbbm{1}_{t-1}&\mathbbm{1}_{t}&\mathbbm{1}_{s-1}&\mathbbm{1}_{t}&\mathbbm{1}_{s}&\mathbbm{1}_{s}&\mathbbm{1}_{t-1}\\\
-1&1&1&1&\mathbbm{1}_{s-1}&\mathbbm{1}_{t-1}&\mathbbm{1}_{t}&\mathbbm{1}_{s-1}&-\mathbbm{1}_{t}&-\mathbbm{1}_{s}&-\mathbbm{1}_{s}&-\mathbbm{1}_{t-1}\\\
-1&-1&1&1&\mathbbm{1}_{s-1}&\mathbbm{1}_{t-1}&-\mathbbm{1}_{t}&-\mathbbm{1}_{s-1}&\mathbbm{1}_{t}&\mathbbm{1}_{s}&-\mathbbm{1}_{s}&-\mathbbm{1}_{t-1}\\\
-1&-1&-1&1&\mathbbm{1}_{s-1}&-\mathbbm{1}_{t-1}&\mathbbm{1}_{t}&-\mathbbm{1}_{s-1}&\mathbbm{1}_{t}&-\mathbbm{1}_{s}&\mathbbm{1}_{s}&-\mathbbm{1}_{t-1}\end{array}$
(11)

or

$\displaystyle\begin{array}[]{rrrrrrrrrrrrrrrrrr}1&1&1&1&\mathbbm{1}_{s}&\mathbbm{1}_{t-1}&\mathbbm{1}_{t-1}&\mathbbm{1}_{s}&\mathbbm{1}_{t-1}&\mathbbm{1}_{s}&\mathbbm{1}_{s}&\mathbbm{1}_{t-1}\\\
-1&1&1&-1&\mathbbm{1}_{s}&\mathbbm{1}_{t-1}&\mathbbm{1}_{t-1}&\mathbbm{1}_{s}&-\mathbbm{1}_{t-1}&-\mathbbm{1}_{s}&-\mathbbm{1}_{s}&-\mathbbm{1}_{t-1}\\\
-1&-1&1&1&\mathbbm{1}_{s}&\mathbbm{1}_{t-1}&-\mathbbm{1}_{t-1}&-\mathbbm{1}_{s}&\mathbbm{1}_{t-1}&\mathbbm{1}_{s}&-\mathbbm{1}_{s}&-\mathbbm{1}_{t-1}\\\
-1&1&-1&1&\mathbbm{1}_{s}&-\mathbbm{1}_{t-1}&\mathbbm{1}_{t-1}&-\mathbbm{1}_{s}&\mathbbm{1}_{t-1}&-\mathbbm{1}_{s}&\mathbbm{1}_{s}&-\mathbbm{1}_{t-1}\end{array}$
(16)

for some $s,t$ with $s+t=\frac{n}{4}$. The quadruples of rows given in (11)
and (16) are of type $\min\\{s,t\\}$. Letting

$U=\left[\begin{array}[]{rrrr}0&0&0&-1\\\ 1&0&0&0\\\ 0&1&0&0\\\ 0&0&1&0\\\
\end{array}\right]$

and
$V=\mathrm{diag}(-I_{s-1},I_{t-1},-I_{t},I_{s-1},-I_{t},I_{s},-I_{s},I_{t-1})$
and multiplying the quadruple of rows given in (11) by $U$ from the left and
by $U^{\top}\oplus V$ from the right, we may assume that $s\geqslant t$. This
means the quadruple of rows given in (11) is of type $t$. We define the skew
type of the quadruples of rows given in (11) as $(t,1)$. Also, we define the
skew type of the quadruples of rows given in (16) as $(s,2)$ if $s\leqslant t$
and as $(t,0)$ if $s>t$. Note that skew types $(0,0)$ and $(0,1)$ are not
admissible. A skew Hadamard matrix is of skew type $(t,\varepsilon)$ if it is
a Hadamard matrix of type $t$ containing a quadruple of rows of skew type
$(t,\varepsilon)$ and no quadruples of rows of skew type
$(t,\varepsilon^{\prime})$ with $\varepsilon^{\prime}<\varepsilon$.

For a quadruple $\\{i,j,k,\ell\\}$ of rows of a skew Hadamard matrix
$H=[h_{ij}]$, define

$Q_{ijk\ell}=\sum_{r\in\\{i,j,k,\ell\\}}h_{ir}h_{jr}h_{kr}h_{\ell r}.$

Note that $Q_{ijk\ell}=0$ if and only if the $4\times 4$ principal submatrix
corresponding to the quadruple $\\{i,j,k,\ell\\}$ is SH-equivalent to $A$,
given in (6). Moreover, the parameter $P_{ijk\ell}Q_{ijk\ell}$ is invariant
under negation and permutation of rows, and simultaneously, negation and
permutation of columns The following proposition describes the skew type of a
quadruple of rows of a skew Hadamard matrix. It is easily proved using the
forms (11) and (16) and noting that $P_{ijk\ell}Q_{ijk\ell}$ is an invariant
parameter.

######  Proposition 2.1.

Let $H$ be a skew Hadamard matrix of order $n$. Then, the following statements
hold for every quadruple $\\{i,j,k,\ell\\}$ of rows of $H$ of type $t$.

  * (i)

The skew type of $\\{i,j,k,\ell\\}$ is $(t,0)$ if and only of
$P_{ijk\ell}Q_{ijk\ell}<0$.

  * (ii)

The skew type of $\\{i,j,k,\ell\\}$ is $(t,1)$ if and only of $Q_{ijk\ell}=0$.

  * (iii)

The skew type of $\\{i,j,k,\ell\\}$ is $(t,2)$ if and only of $Q_{ijk\ell}\neq
0$ and $P_{ijk\ell}Q_{ijk\ell}\geqslant 0$.

The following lemma shows the relation between the skew type of a quadruple of
rows in a skew Hadamard matrix $H$ and the skew type of their corresponding
rows in $H^{\top}$.

######  Lemma 2.2.

Let $H$ be a skew Hadamard matrix of order $n$. Then, the following statements
hold.

  * (i)

If the quadruple $\\{i,j,k,\ell\\}$ of rows of $H$ is of skew type $(t,0)$,
then the quadruple $\\{i,j,k,\ell\\}$ of rows of $H^{\top}$ is of skew type
$(t-1,2)$.

  * (ii)

If the quadruple $\\{i,j,k,\ell\\}$ of rows of $H$ is of skew type $(t,1)$,
then the quadruple $\\{i,j,k,\ell\\}$ of rows of $H^{\top}$ is also of skew
type $(t,1)$.

  * (iii)

If the quadruple $\\{i,j,k,\ell\\}$ of rows of $H$ is of skew type $(s,2)$,
then the quadruple $\\{i,j,k,\ell\\}$ of rows of $H^{\top}$ is of skew type

$\left\\{\begin{array}[]{ll}(s-1,2)&\mbox{ if }s=\frac{n}{8}\mbox{,}\\\
\vspace{-1mm}\\\ (s,2)&\mbox{ if }s=\frac{n-4}{8}\mbox{,}\\\ \vspace{-1mm}\\\
(s+1,2)&\mbox{ if }s=\frac{n}{8}-1\mbox{,}\\\ \vspace{-1mm}\\\ (s+1,0)&\mbox{
}\mbox{otherwise.}\end{array}\right.$

###### Proof.

To prove (i), let $\\{i,j,k,\ell\\}$ be a quadruple of rows of $H$ of skew
type $(t,0)$ and let $s=\frac{n}{4}-t$. By (16), we have
$T_{ijk\ell}=\frac{n-|P_{ijk\ell}|}{8}=\frac{4(s+t)-4|s-t+2|}{8}=t-1$ in
$H^{\top}$, by noting that $s>t$. Since $P_{ijk\ell}Q_{ijk\ell}>0$ in
$H^{\top}$, it follows from Proposition 2.1(iii) that the quadruple
$\\{i,j,k,\ell\\}$ of rows of $H^{\top}$ is of skew type $(t-1,2)$.

For (ii), let $\\{i,j,k,\ell\\}$ be a quadruple of rows of $H$ of skew type
$(t,1)$ and let $s=\frac{n}{4}-t$. By (11), we have
$T_{ijk\ell}=\frac{n-|P_{ijk\ell}|}{8}=\frac{4(s+t)-4|s-t|}{8}=t$ in
$H^{\top}$, by noting that $s\geqslant t$. Since $Q_{ijk\ell}=0$ in
$H^{\top}$, it follows from Proposition 2.1(ii) that the quadruple
$\\{i,j,k,\ell\\}$ of rows of $H^{\top}$ is of skew type $(t,1)$.

Finally, to establish (iii), let $\\{i,j,k,\ell\\}$ be a quadruple of rows of
$H$ of skew type $(s,2)$ and let $t=\frac{n}{4}-s$. As $s\leqslant t$, it
follows from (16) that

$T_{ijk\ell}=\frac{n-|P_{ijk\ell}|}{8}=\frac{4(s+t)-4|s-t+2|}{8}=\left\\{\begin{array}[]{ll}s-1&\mbox{
if }s=\frac{n}{8}\mbox{,}\\\ \vspace{-1mm}\\\ s&\mbox{ if
}s=\frac{n-4}{8}\mbox{,}\\\ \vspace{-1mm}\\\ s+1&\mbox{
}\mbox{otherwise,}\end{array}\right.$

in $H^{\top}$. If $s\in\\{\frac{n}{8},\frac{n-4}{8},\frac{n}{8}-1\\}$, then
$P_{ijk\ell}Q_{ijk\ell}\geqslant 0$ in $H^{\top}$, and otherwise,
$P_{ijk\ell}Q_{ijk\ell}<0$ in $H^{\top}$. As $Q_{ijk\ell}\neq 0$ in
$H^{\top}$, the result follows from Proposition 2.1. ∎

We recall the following two known lemmas on the types of a Hadamard matrix and
its transpose.

######  Lemma 2.3 ([12, 14, 19]).

There is no Hadamard matrix of order $n$ and type $\lfloor\frac{n}{8}\rfloor$
if $n\notin\\{4,12\\}$. There is no Hadamard matrix of order $n$ and type $0$
if $n\equiv 4\pmod{8}$.

######  Lemma 2.4 ([12]).

Let $H$ be a Hadamard matrix of order $n$ having a quadruple of rows of type
$1$. Then, the type of $H^{\top}$ is $0$ if $n\equiv 0\pmod{8}$ and $1$ if
$n\equiv 4\pmod{8}$.

######  Theorem 2.5.

The following statements hold.

  * (i)

There is no skew Hadamard matrix of skew type $(2,0)$.

  * (ii)

There is no skew Hadamard matrix of order $n$ and skew type $(1,0)$ if
$n\equiv 4\pmod{8}$.

  * (iii)

There is no skew Hadamard matrix of order $n$ and skew type $(1,1)$ or $(1,2)$
if $n\equiv 0\pmod{8}$.

###### Proof.

Let $H$ be a skew Hadamard matrix of order $n$. Suppose the skew type of $H$
is $(2,0)$. Then, Lemma 2.2(i) implies that $H^{\top}$ has a quadruple of rows
of skew type $(1,2)$ and so $H^{\top}$ is of type $0$ or $1$. If $H^{\top}$ is
of type $1$, then, by applying Lemma 2.4 for $H^{\top}$, we deduce that $H$ is
of type $0$ or $1$ which is a contradiction, since the skew type of $H$ is
$(2,0)$. Therefore, $H^{\top}$ is of type $0$ and hence $n\equiv 0\pmod{8}$ by
Lemma 2.3. Since the skew type of $H^{\top}$ is $(0,2)$, Lemma 2.2(iii)
implies that $H$ has a quadruple of rows of type $1$, which is again a
contradiction, as the skew type of $H$ is $(2,0)$.

Now, let $n\equiv 4\pmod{8}$. If $H$ is of skew type $(1,0)$, then Lemma
2.2(i) implies that $H^{\top}$ has a quadruple of rows of skew type $(0,2)$
and so $H^{\top}$ is of type $0$, contradicting Lemma 2.3.

Finally, suppose that $n\equiv 0\pmod{8}$ and $H$ is of skew type $(1,1)$ or
$(1,2)$. It follows from Lemma 2.3 that $n\neq 8$. Also, Lemma 2.4 yields that
the type of $H^{\top}$ is $0$, so $H^{\top}$ has a quadruple of skew type
$(0,2)$ rows. So, it follows from Lemma 2.2(iii) that $H$ has a quadruple of
rows of skew type $(1,0)$. As $H$ is of type $1$, we conclude that the skew
type of $H$ is $(1,0)$, a contradiction. ∎

To proceed, we must recall the following result on the type of quadruples of
rows of a Hadamard matrix.

######  Lemma 2.6 ([19]).

Let $H$ be a Hadamard matrix of order $4m$. Fix three rows of $H$ and let
$\kappa_{t}$ be the number of other rows of type $t$ with these three rows.
Then,

$\sum_{t=0}^{\left\lfloor\frac{m}{2}\right\rfloor}\kappa_{t}(m-2t)^{2}=m^{2}.$

######  Lemma 2.7.

Let $H$ be a skew Hadamard matrix of order $n$. Then, the number of all
$4\times 4$ principal submatrices of $H$ which are SH-equivalent to $A$ is
$n(n-1)(n-2)(n-4)/32$, where $A$ is introduced in (6).

###### Proof.

For simplicity, let $n=4m$. It is easy to see that, by a sequence of row
negations and row permutations, and simultaneously, column negations and
column permutations, every triple of rows of $H$ may be transformed to the
form

$\displaystyle\begin{array}[]{rrrrrrrrrrrrrrrrrr}1&1&1&\mathbbm{1}_{m-1}&\mathbbm{1}_{m-1}&\mathbbm{1}_{m-1}&\mathbbm{1}_{m}\\\
-1&1&1&\mathbbm{1}_{m-1}&\mathbbm{1}_{m-1}&-\mathbbm{1}_{m-1}&-\mathbbm{1}_{m}\\\
-1&-1&1&\mathbbm{1}_{m-1}&-\mathbbm{1}_{m-1}&-\mathbbm{1}_{m-1}&\mathbbm{1}_{m}\end{array}.$

One may straightforwardly check that $Q_{123\ell}=0$ for the quadruple of rows
$\\{1,2,3,\ell\\}$ if and only if $\ell\in\\{4,\ldots,3m\\}$. This means that
the quadruple of rows $\\{1,2,3,\ell\\}$ forms a $4\times 4$ principal
submatrix of $H$ which is SH-equivalent to $A$ if and only if
$\ell\in\\{4,\ldots,3m\\}$. Therefore, every triple of rows of $H$ is
contained in exactly $3m-3$ quadruples of rows whose corresponding $4\times 4$
principal submatrix is SH-equivalent to $A$. By double counting, we derive
that the number of all quadruples of rows of $H$ whose corresponding $4\times
4$ principal submatrices are SH-equivalent to $A$ is equal to
$(3m-3){{n}\choose{3}}/4=n(n-1)(n-2)(n-4)/32$, as desired. ∎

######  Theorem 2.8.

There is no skew Hadamard matrix of order $36$ and skew type $(3,2)$.

###### Proof.

Suppose by way of contradiction that $H$ is a skew Hadamard matrix of order
$36$ and skew type $(3,2)$. So, every quadruple of rows of $H$ has type $3$ or
$4$. Using the notation of Lemma 2.6, we have $9\kappa_{3}^{2}+\kappa_{4}=81$.
Since $\kappa_{3}+\kappa_{4}=33$, we obtain that $\kappa_{3}=6$ and
$\kappa_{4}=27$. By double counting, we derive that the number of all
quadruples of rows of $H$ of type $3$ is $6{{36}\choose{3}}/4=10710$ and
therefore, the number of all quadruples of rows of $H$ of type $4$ is
${{36}\choose{4}}-10710=48195$.

The possible skew types for any quadruple of rows of $H$ are $(3,2)$, $(4,0)$,
$(4,1)$, $(4,2)$. Using Lemma 2.2, the possible skew types for any quadruple
of rows of $H^{\top}$ are $(3,2)$, $(4,0)$, $(4,1)$, $(4,2)$ and therefore,
$H^{\top}$ is of skew type $(3,2)$. Since the number of all quadruples of rows
of $H$ of type $3$ is $10710$, the number of all quadruples of rows of $H$ of
skew type $(3,2)$ is $10710$. Also, Lemma 2.2 implies that the number of all
quadruples of rows of $H$ of skew type $(4,0)$ is equal to the number of all
quadruples of rows of $H^{\top}$ of skew type $(3,2)$. Therefore, the number
of all quadruples of rows of $H$ of skew type $(4,0)$ is $10710$.

The number of all quadruples of rows of $H$ of skew type $(4,1)$ is equal to
the number of all quadruples of rows of $H$ whose corresponding $4\times 4$
principal submatrices are SH-equivalent to $A$ which equals $42840$, by Lemma
2.7. Now, the number of all quadruples of rows of $H$ of skew type $(4,2)$ is
$48195-(10710+42840)<0$, a contradiction. ∎

We end this section with the following lemma, which is crucial in our
computational search.

######  Lemma 2.9.

Let $H$ be a Hadamard matrix of order $n$ whose first four rows are in the
form (5). Let $(h_{1},\ldots,h_{n})$ be a row of $H$ other than the first four
rows and define

$x_{1}=\sum_{i=1}^{s}h_{i},\qquad x_{2}=\sum_{i=s+2t+1}^{2s+2t}h_{i},\qquad
x_{3}=\sum_{i=2s+3t+1}^{3s+3t}h_{i},\qquad
x_{4}=\sum_{i=3s+3t+1}^{4s+3t}h_{i}$

and

$y_{1}=\sum_{i=s+1}^{s+t}h_{i},\qquad y_{2}=\sum_{i=s+t+1}^{s+2t}h_{i},\qquad
y_{3}=\sum_{i=2s+2t+1}^{2s+3t}h_{i},\qquad y_{4}=\sum_{i=4s+3t+1}^{n}h_{i}.$

Then, $x_{1}+x_{2}+x_{3}-x_{4}=2s-\frac{n}{2}\pmod{4}$ and
$y_{1}+y_{2}+y_{3}-y_{4}=2t-\frac{n}{2}\pmod{4}$.

###### Proof.

Since $(h_{1},\ldots,h_{n})$ is orthogonal to the four rows given in (5), we
obtain that

$\displaystyle x_{1}+y_{1}+y_{2}+x_{2}+y_{3}+x_{3}+x_{4}+y_{4}$
$\displaystyle=0,$ (18) $\displaystyle
x_{1}+y_{1}+y_{2}+x_{2}-y_{3}-x_{3}-x_{4}-y_{4}$ $\displaystyle=0,$ (19)
$\displaystyle x_{1}+y_{1}-y_{2}-x_{2}+y_{3}+x_{3}-x_{4}-y_{4}$
$\displaystyle=0,$ (20) $\displaystyle
x_{1}-y_{1}+y_{2}-x_{2}+y_{3}-x_{3}+x_{4}-y_{4}$ $\displaystyle=0.$ (21)

By subtracting (21) from the equality obtained by taking the sum of (18)–(20),
we conclude that $x_{1}+x_{2}+x_{3}-x_{4}=-2y_{1}$. Since $2y_{1}=2t\pmod{4}$
and $s+t=\frac{n}{4}$, one deduces that
$x_{1}+x_{2}+x_{3}-x_{4}=2s-\frac{n}{2}\pmod{4}$. Also, by summing (18)–(21),
we get $y_{1}+y_{2}+y_{3}-y_{4}=-2x_{1}$. Since $2x_{1}=2s\pmod{4}$ and
$s+t=\frac{n}{4}$, we get $y_{1}+y_{2}+y_{3}-y_{4}=2t-\frac{n}{2}\pmod{4}$,
the result follows. ∎

##  3 Search using skew types

In this section, we describe our search method for skew Hadamard matrices,
which uses the notion of skew types. Let $H$ be a skew Hadamard matrix of
order $n$ and skew type $(t,\varepsilon)$. By Lemma 2.2, one can easily see
that if $\varepsilon=0$, then $H^{\top}$ is of skew type $(t-1,2)$. Therefore,
the classification of skew Hadamard matrices of skew type $(t,0)$ can be
extracted from the classification of skew Hadamard matrices of skew type
$(t-1,2)$. Hence, we may assume that $\varepsilon=1,2$. Using the forms (11),
(16) and by a permutation on columns, we may assume that the first four rows
of $H$ are of the form

$\displaystyle\begin{array}[]{rrrrrrrrrrrrrrrrrr}1&1&1&1&\mathbbm{1}_{t-1}&\mathbbm{1}_{t}&\mathbbm{1}_{t}&\mathbbm{1}_{t-1}&\mathbbm{1}_{s-1}&\mathbbm{1}_{s-1}&\mathbbm{1}_{s}&\mathbbm{1}_{s}\\\
-1&1&1&1&\mathbbm{1}_{t-1}&\mathbbm{1}_{t}&-\mathbbm{1}_{t}&-\mathbbm{1}_{t-1}&\mathbbm{1}_{s-1}&\mathbbm{1}_{s-1}&-\mathbbm{1}_{s}&-\mathbbm{1}_{s}\\\
-1&-1&1&1&\mathbbm{1}_{t-1}&-\mathbbm{1}_{t}&\mathbbm{1}_{t}&-\mathbbm{1}_{t-1}&\mathbbm{1}_{s-1}&-\mathbbm{1}_{s-1}&\mathbbm{1}_{s}&-\mathbbm{1}_{s}\\\
-1&-1&-1&1&-\mathbbm{1}_{t-1}&\mathbbm{1}_{t}&\mathbbm{1}_{t}&-\mathbbm{1}_{t-1}&\mathbbm{1}_{s-1}&-\mathbbm{1}_{s-1}&-\mathbbm{1}_{s}&\mathbbm{1}_{s}\end{array}$
(26)

or

$\displaystyle\begin{array}[]{rrrrrrrrrrrrrrrrrr}1&1&1&1&\mathbbm{1}_{t}&\mathbbm{1}_{t}&\mathbbm{1}_{t}&\mathbbm{1}_{t}&\mathbbm{1}_{s-1}&\mathbbm{1}_{s-1}&\mathbbm{1}_{s-1}&\mathbbm{1}_{s-1}\\\
-1&1&1&-1&\mathbbm{1}_{t}&\mathbbm{1}_{t}&-\mathbbm{1}_{t}&-\mathbbm{1}_{t}&\mathbbm{1}_{s-1}&\mathbbm{1}_{s-1}&-\mathbbm{1}_{s-1}&-\mathbbm{1}_{s-1}\\\
-1&-1&1&1&\mathbbm{1}_{t}&-\mathbbm{1}_{t}&\mathbbm{1}_{t}&-\mathbbm{1}_{t}&\mathbbm{1}_{s-1}&-\mathbbm{1}_{s-1}&\mathbbm{1}_{s-1}&-\mathbbm{1}_{s-1}\\\
-1&1&-1&1&\mathbbm{1}_{t}&-\mathbbm{1}_{t}&-\mathbbm{1}_{t}&\mathbbm{1}_{t}&-\mathbbm{1}_{s-1}&\mathbbm{1}_{s-1}&\mathbbm{1}_{s-1}&-\mathbbm{1}_{s-1}\end{array}$
(31)

where $s=\frac{n}{4}-t\geqslant t$.

By applying Lemma 2.9, if the first four rows of $H$ are of the form (26),
then the column $4t+2$ of $H$ is uniquely determined from the columns
$1,3,5,\dots,4t+1$ of $H$ and hence the row $4t+2$ of $H$ is uniquely
determined from the rows $1,3,5,\dots,4t+1$ of $H$. Also, using Lemma 2.9, if
the first four rows of $H$ are of the form (31), then the column $4t+4$ of $H$
is uniquely determined from the columns $5,\dots,4t+3$ of $H$ and thus the row
$4t+4$ of $H$ is uniquely determined from the rows $5,\dots,4t+3$ of $H$.

To classify skew Hadamard matrices of skew type $(t,\varepsilon)$ for
$\varepsilon=1,2$, we first fix the first four rows according to one of the
forms (26) and (31). We then use a backtrack algorithm to construct the next
rows one by one. By the above paragraph, when we reach the row $4t+2$ or
$4t+4$, depending on the form of the first four rows, we have pruning criteria
for the backtracking. We compute this row from the preceding rows. If there is
no solution, we backtrack. Otherwise, we proceed to the next step.

During the backtracking algorithm, we perform some checks. Regarding the
number of solutions permitted in the primary steps of backtracking, we do the
SH-equivalence test on the obtained partial skew Hadamard matrices. Also, we
check for the skew type of partial skew Hadamard matrices whenever it is
computationally not time-consuming.

After we obtain all solutions, we need to check them for the SH-equivalence.
It could take a long time because of the large number of solutions. To do
this, we partition the set of solutions according to their skew profiles,
which count the number of quadruples of rows of different skew types. Then, in
each class, which is usually much smaller, we do the SH-equivalence test.

In order $36$, in view of Lemma 2.3, there is no skew Hadamard matrix of skew
type $(0,\varepsilon)$ and $(4,\varepsilon)$ for any $\epsilon$. By Lemma 2.5,
there is no such matrix of skew type $(1,0)$ and $(2,0)$. We also have no
matrix of type $(3,2)$, using Lemma 2.8. So, we need to find solutions for
skew types $(1,1)$, $(1,2)$, $(2,1)$, $(2,2)$, $(3,0)$ and $(3,1)$. We have
run the program for skew types $(1,1)$, $(1,2)$ and $(2,1)$. The two skew
types $(1,1)$ and $(1,2)$ are fast enough to be run on a single PC for a few
days. The computation time for the skew type $(2,1)$ was about twenty days on
a cluster of 48 cores. We also implemented some parts of the algorithm twice
with different codes. The remaining cases $(2,2)$, $(3,0)$ and $(3,1)$ need
much more computational resources. We also run the program for all orders
smaller than $36$. The results, which are given in Table 2, confirm the
numbers shown in Table 1. The obtained list of skew Hadamard matrices of order
$36$ is available electronically at [9]. We summarize our main result in the
following proposition.

######  Proposition 3.1.

There are at least $157132$ SH-inequivalent skew Hadamard matrices of order
$36$.

Order Skew Type | (0, 2) | (1, 0) | (1, 1) | (1, 2) | (2, 1) | (2, 2) | (3, 0) | (3, 1)
---|---|---|---|---|---|---|---|---
4 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0
8 | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 0
12 | 0 | 0 | 1 | 0 | 0 | 0 | 0 | 0
16 | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 0
20 | 0 | 0 | 1 | 1 | 0 | 0 | 0 | 0
24 | 14 | 1 | 0 | 0 | 1 | 0 | 0 | 0
28 | 0 | 0 | 43 | 21 | 0 | 1 | 0 | 0
32 | 6903 | 283 | 0 | 0 | 40 | 0 | 0 | 1
36 | 0 | 0 | 23260 | 123326 | 10546 | ? | ? | ?
Table 2: The number of SH-inequivalent skew Hadamard matrices of orders up to
$36$ for all possible skew types.

##  4 Search through the canonical method

This section uses the orderly generation method to classify skew Hadamard
matrices of order $36$. The technique was independently introduced in [7] and
[22]. The method consists of two parallel subroutines. These are the
construction of objects and the rejection of equivalence copies. The most
natural and widely used method for the construction phase is backtracking,
which has quite an old history. We use the notion of ‘canonical form’ to
reject equivalence copies. A canonical form is a special representative for
each equivalence class. The canonical forms are constructed step by step using
backtracking, and the canonicity test is carried out at each step.

Here, we present the canonical form for skew Hadamard matrices. For a skew
Hadamard matrix $H=[h_{ij}]$ of order $n$, assign a vector

$v(H)=\big{(}h_{1,2},h_{1,3},h_{2,3},h_{1,4},h_{2,4},h_{3,4},\ldots,h_{n-1,n}\big{)}.$

We say that $H$ is in the canonical form if $v(P^{\top}HP)\preccurlyeq v(H)$
for any signed permutation matrix $P$, where $\preccurlyeq$ denotes the
lexicographical order.

We used an orderly algorithm with backtracking and the above canonical form to
generate SH-equivalence classes of skew Hadamard matrices of order $36$. Using
this orderly algorithm, we first find all solutions for the $14\times 14$
leading principal submatrices in canonical form. In total, there are
$80122802$ solutions for $14\times 14$ matrices. For each solution, we find
all $14\times 36$ partial Hadamard matrices constituting the first $14$ rows
of desired skew Hadamard matrix of order $36$ using an exhaustive search.
Then, for each obtained $14\times 36$ partial Hadamard matrix and
$i=15,\ldots,36$, we complete every row $i$ by solving a system of linear
equations coming from the orthogonality row $i$ to the first 14 rows. To speed
up the computation, we solve the system of linear equations in modulo $2$.
Finally, for the complete solutions, we do the canonicity test. We have not
been able to run the computation for all $80122802$ solutions. It is done for
about $10$ million of them. The approximate computation time was three months
on a cluster of 160 cores. We found $157132$ skew Hadamard matrices of order
$36$ up to the SH-equivalence. The list matches the list obtained in the
previous section.

Denote by $J_{k}$ the skew-symmetric matrix of order $k$ whose all entries on
and above the main diagonal are $1$. For example, the matrix $A$ given in (6)
equals $J_{4}$. The number of solutions for the $8\times 8$ leading principal
submatrices of skew Hadamard matrices of order $36$ in the canonical form is
$9$. Our computation shows that only one of these $9$ solutions, $J_{8}$, is
extendable to skew Hadamard matrices. So, every skew Hadamard matrix of order
$36$ in canonical form contains $J_{8}$. Note that $J_{8}$ has $11$ feasible
extensions to $9\times 9$ submatrices of skew Hadamard matrices of order $36$.
Let us show these solutions by $S_{0}=J_{9}$, $S_{1},\ldots,S_{10}$ from the
largest to smallest in the canonical form. Our results show that
$S_{6},\ldots,S_{10}$ have no extensions. For $S_{0}$ and $S_{1}$, we found
the complete list of extensions. However, the search is not complete for
$S_{2}$. Also, we have not done anything for $S_{3}$, $S_{4}$ and $S_{5}$.

Inspired by the previous paragraph, we may ask for the maximum number $k(n)$
such that $J_{k(n)}$ appears as a principal submatrix of a skew Hadamard
matrix of order $n$. This is related to a problem posed in [6]: What is the
maximum number $t(n)$ such that every tournament of order $n$ contains the
transitive tournament of order $t(n)$ as a subtournament? The following
theorem presents an upper bound on $k(n)$, which is especially tight for
$n=36$.

######  Theorem 4.1.

For any $n$,

$k(n)\leqslant\left\lfloor\frac{\pi}{2\cot^{-1}\big{(}\sqrt{n-1}\big{)}}\right\rfloor.$

###### Proof.

Let $H$ be a skew Hadamard matrix of order $n$. It is easy to see that
$A=\mathbbm{i}(I-H)$ is a Hermitian matrix with $\frac{n}{2}$ eigenvalues
$\sqrt{n-1}$ and $\frac{n}{2}$ eigenvalues $-\sqrt{n-1}$, where $\mathbbm{i}$
is the unit imaginary complex number. Let $k=k(n)$ and
$B=\mathbbm{i}(I-J_{k})$. As $J_{k}$ is a principal submatrix of $H$, $B$ is a
principal submatrix of $A$. Since $B$ is a negacirculant matrix, it
straightforwardly follows from a result of [15] that the eigenvalues of $B$
are $\cot(\frac{\pi\ell}{2k})$ for $\ell=1,3,\ldots,2k-1$. By the interlacing
theorem, the largest eigenvalue of $B$ is less than or equal to the largest
eigenvalue of $A$, that is, $\cot(\frac{\pi}{2k})\leqslant\sqrt{n-1}$. This
yields the desired inequality. ∎

######  Remark 4.2.

Whether or not our list of skew Hadamard matrices of order $36$ is complete is
a question. Referring to Table 2, we strongly feel that no skew Hadamard
matrices of order $36$ exist with skew types $(2,2)$, $(3,0)$ and $(3,1)$.

##  5 Ternary codes

Self-dual codes are one of the most interesting classes of codes. This
interest is justified by many combinatorial objects and algebraic objects
related to self-dual codes, see, for example, [21]. In this section, we
classify ternary self-dual codes constructed from the $157132$ skew Hadamard
matrices of order $36$ given in Table 2. We also study unimodular lattices and
$1$-designs built from the ternary near-extremal self-dual codes.

###  5.1 Ternary self-dual codes

Let $\mathbbmsl{F}_{3}=\\{0,1,2\\}$ denote the finite field of order $3$. A
ternary $[n,k]$ code $C$ is a $k$-dimensional vector subspace of
$\mathbbmsl{F}_{3}^{n}$. The parameters $n$ and $k$ are called the length and
dimension of $C$, respectively. A ternary self-dual code $C$ of length $n$ is
a ternary $[n,n/2]$ code satisfying $C=C^{\perp}$, where $C^{\perp}$ is the
dual code of $C$ under the standard inner product. It was shown in [16] that
the minimum weight $d$ of a ternary self-dual code of length $n$ is bounded by
$d\leqslant 3\lfloor n/12\rfloor+3$. If $d=3\lfloor n/12\rfloor+3$
(respectively, $d=3\lfloor n/12\rfloor$), then the code is called extremal
(respectively, near-extremal). For length $36$, the Pless symmetry code is a
currently known extremal ternary self-dual code.

Two ternary codes $C$ and $C^{\prime}$ are said to be equivalent if there
exists a monomial matrix $P$ over $\mathbbmsl{F}_{3}$ with
$C^{\prime}=\\{xP\,|\,x\in C\\}$, otherwise, they are called inequivalent.
Here, a monomial matrix is a matrix with entries in $\mathbbmsl{F}_{3}$ with
exactly one nonzero entry in each row and each column. Let $H$ be a Hadamard
matrix of order $n$. Throughout this section, let $C(H)$ denote the ternary
code generated by the rows of $H$, where the entries of the matrix are
regarded as elements of $\mathbbmsl{F}_{3}$. It is trivial that $C(H)$ and
$C(K)$ are equivalent if $H$ and $K$ are SH-equivalent skew Hadamard matrices.

Although the proof of the following lemma is somewhat trivial, we give it for
the sake of completeness.

######  Lemma 5.1.

Let $H$ be a skew Hadamard matrix of order $n$. If $n\equiv 0\pmod{12}$, then
$C(H)$ is self-dual.

###### Proof.

Since $n\equiv 0\pmod{12}$, $HH^{\top}\equiv 0_{n}\pmod{3}$, where $0_{n}$
denotes the $n\times n$ zero matrix. This implies that $C(H)\subseteq
C(H)^{\perp}$. From $\dim C(H)+\dim C(H)^{\perp}=n$, we find that $\dim
C(H)\leqslant\frac{n}{2}$. On the other hand, since $H+H^{\top}=2I$,

$n=\operatorname{rank}_{3}(2I)=\operatorname{rank}_{3}\left(H+H^{\top}\right)\leqslant\operatorname{rank}_{3}(H)+\operatorname{rank}_{3}\left(H^{\top}\right)=2\operatorname{rank}_{3}(H),$

where $\operatorname{rank}_{3}(X)$ denotes the rank of the matrix $X$ regarded
over $\mathbbmsl{F}_{3}$. Thus,
$\operatorname{rank}_{3}(H)\geqslant\frac{n}{2}$ which yields that $\dim
C(H)\geqslant\frac{n}{2}$. Therefore, $\dim C(H)=\frac{n}{2}$ and so
$C(H)=C(H)^{\perp}$, meaning that $C(H)$ is self-dual. ∎

We describe how to classify ternary self-dual codes of length $36$ constructed
from the $157132$ skew Hadamard matrices of order $36$ given in Table 2. By
comparing pairs $(A_{6},A_{9})$, the $157132$ ternary self-dual codes are
divided into $829$ classes, where $A_{w}$ denotes the number of codewords of
weight $w$. For each class, to test the equivalence of ternary self-dual
codes, we used the method given in Section 7.3.3 of [10] as follows. For a
ternary self-dual code $C$, define the vertex-colored graph
$\mathnormal{\Gamma}(C)$ with vertex set
$C_{9}\cup(\\{1,2,\dots,36\\}\times(\mathbbmsl{F}_{3}\setminus\\{0\\}))$ and
edge set $\\{\\{c,(i,c_{i})\\}\,|\,c=(c_{1},\ldots,c_{36})\in
C_{9},i\in\\{1,\ldots,36\\}\text{ and }c_{i}\neq
0\\}\cup\\{\\{(i,y),(i,2y)\\}\,|\,i\in\\{1,\ldots,36\\}\text{ and
}y\in\mathbbmsl{F}_{3}\setminus\\{0\\}\\}$, where $C_{9}$ denotes the set of
codewords of weight $9$ in $C$. Suppose that two ternary self-dual codes $C$
and $C^{\prime}$ are generated by codewords of weight $9$. Then, $C$ and
$C^{\prime}$ are equivalent if and only if $\mathnormal{\Gamma}(C)$ and
$\mathnormal{\Gamma}(C^{\prime})$ are isomorphic. We verified that the
$157129$ ternary self-dual codes of length $36$ are generated by codewords of
weight $9$ and the remaining three codes are not generated by codewords of
weight $9$. For the three codes, we have that $(A_{6},A_{9})=(36,464)$,
$(96,224)$ and $(96,1520)$. Thus, each of the three codes is inequivalent to
the other codes. We used nauty [17] for vertex-colored graph isomorphism
testing corresponding to the $157129$ ternary self-dual codes. Then, we have
the following classification of ternary self-dual codes of length $36$
constructed from the $157132$ skew Hadamard matrices of order $36$.

######  Proposition 5.2.

$153979$ inequivalent ternary self-dual codes of length $36$ are constructed
from the $157132$ skew Hadamard matrices of order $36$ given in Table 2.
$20848$ of them have minimum weight $9$ and the others have minimum weight
$6$.

The inequivalence of the above $153979$ ternary self-dual codes was verified
independently by Magma [2]. The obtained list of the ternary self-dual codes
of length $36$ is available electronically at [9].

###  5.2 Ternary near-extremal self-dual codes

Here, we concentrate on the $20848$ ternary near-extremal self-dual codes
given in Proposition 5.2. Recently, a restriction on the weight enumerators of
ternary near-extremal self-dual codes of lengths divisible by $12$ has been
given in [1]. The possible weight enumerators of ternary near-extremal self-
dual codes of length $36$ are

$\displaystyle W_{\alpha}$ $\displaystyle=1+\alpha
y^{9}+(42840-9\alpha)y^{12}+(1400256+36\alpha)y^{15}+(18452280-84\alpha)y^{18}$
$\displaystyle+(90370368+126\alpha)y^{21}+(162663480-126\alpha)y^{24}+(97808480+84\alpha)y^{27}$
$\displaystyle+(16210656-36\alpha)y^{30}+(471240+9\alpha)y^{33}+(888-\alpha)y^{36},$

where $\alpha$ is an integer with $1\leqslant\alpha\leqslant 888$ and
$\alpha\equiv 0\pmod{8}$. For the ternary near-extremal self-dual codes given
in Proposition 5.2, $\alpha$ in their weight enumerators $W_{\alpha}$ are
listed in Table 3. In addition, the numbers $N_{\alpha}$ of the ternary near-
extremal self-dual codes having weight enumerator $W_{\alpha}$ are listed in
Table 3.

$(\alpha,N_{\alpha})$ | $(\alpha,N_{\alpha})$ | $(\alpha,N_{\alpha})$ | $(\alpha,N_{\alpha})$ | $(\alpha,N_{\alpha})$ | $(\alpha,N_{\alpha})$
---|---|---|---|---|---
$(72,1)$ | $(224,25)$ | $(288,1105)$ | $(352,1352)$ | $(416,83)$ | $(480,2)$
$(144,1)$ | $(232,41)$ | $(296,1451)$ | $(360,1123)$ | $(424,49)$ | $(488,2)$
$(168,3)$ | $(240,105)$ | $(304,1580)$ | $(368,870)$ | $(432,34)$ | $(496,1)$
$(176,1)$ | $(248,153)$ | $(312,1860)$ | $(376,679)$ | $(440,26)$ | $(512,1)$
$(184,1)$ | $(256,263)$ | $(320,1750)$ | $(384,466)$ | $(448,9)$ | $(544,1)$
$(200,3)$ | $(264,409)$ | $(328,1894)$ | $(392,342)$ | $(456,15)$ | $(600,1)$
$(208,9)$ | $(272,577)$ | $(336,1762)$ | $(400,200)$ | $(464,8)$ |
$(216,13)$ | $(280,864)$ | $(344,1566)$ | $(408,142)$ | $(472,5)$ |
Table 3: The number $N_{\alpha}$ of the ternary near-extremal self-dual codes.

A (Euclidean) lattice $L\subseteq\mathbbmsl{R}^{n}$ in dimension $n$ is called
unimodular if $L=L^{*}$, where $L^{*}$ is the dual lattice of $L$ under the
standard inner product. The minimum norm of a unimodular lattice $L$ is the
smallest the norm among all nonzero vectors of $L$. The kissing number of $L$
is the number of vectors of minimum norm in $L$. Two lattices $L$ and
$L^{\prime}$ are said to be isomorphic if there exists an orthogonal matrix
$A$ with $L^{\prime}=\\{xA\,|\,x\in L\\}$, otherwise, they are called non-
isomorphic. Let $C$ be a ternary self-dual code of length $n$. It is known
that

$A_{3}(C)=\left\\{\left.\frac{1}{\sqrt{3}}\big{(}x_{1},\ldots,x_{n}\big{)}\in\mathbbmsl{Z}^{n}\,\right|\,\Big{(}x_{1}\pmod{3},\ldots,x_{n}\pmod{3}\Big{)}\in
C\right\\}$

is an unimodular lattice in dimension $n$. This construction of lattices is
well known as Construction A. Let $C$ be a ternary near-extremal self-dual
code of length $36$ having weight enumerator $W_{\alpha}$. Then, it follows
that $A_{3}(C)$ has minimum norm $3$ and kissing number $\alpha+72$. We
verified by Magma [2] that the $20848$ unimodular lattices, which are
constructed from the $20848$ inequivalent ternary near-extremal self-dual
codes given in Proposition 5.2, are non-isomorphic. This also verifies the
inequivalence of the $20848$ ternary near-extremal self-dual codes.

Let $C$ be a ternary near-extremal self-dual code of length $36$ having weight
enumerator $W_{\alpha}$. Then, the supports of codewords of weight $9$ in $C$
form a $1$-$(36,9,\alpha/8)$ design [18]. Let $D_{i}=(X_{i},\mathcal{B}_{i})$
be a $1$-design, where $X_{i}$ is the set of points and $\mathcal{B}_{i}$ is
the collection of blocks of $D_{i}$ for $i=1,2$. Two $1$-designs $D_{1}$ and
$D_{2}$ are said to be isomorphic if there exists a bijection $\phi:X_{1}\to
X_{2}$ such that if we rename every point $x$ of $X_{1}$ by $\phi(x)$ then
$\mathcal{B}_{1}$ is transformed into $\mathcal{B}_{2}$, otherwise, they are
called non-isomorphic. We verified by Magma [2] that these $20848$
$1$-designs, which are constructed from the $20848$ inequivalent ternary near-
extremal self-dual codes given in Proposition 5.2, are non-isomorphic. This
also verifies the inequivalence of the $20848$ ternary near-extremal self-dual
codes.

##  6 Association schemes

In this section, we construct association schemes of order $35$ and class $2$
from skew Hadamard matrices of order $36$. We define association schemes in
matrix form. An association scheme of order $n$ and class $2$ is a set of
nonzero $n\times n$ matrices $\\{A_{0},A_{1},A_{2}\\}$ with entries in
$\\{0,1\\}$ satisfying

  * (i)

$A_{0}=I$;

  * (ii)

$A_{0}+A_{1}+A_{2}=J$, where $J$ denotes the all one matrix;

  * (iii)

For every $i,j\in\\{0,1,2\\}$, $A_{i}^{\top}\in\\{A_{0},A_{1},A_{2}\\}$ and
$A_{i}A_{j}$ is a linear combination of $A_{0},A_{1},A_{2}$.

An association scheme $\mathscr{A}=\\{A_{0},A_{1},A_{2}\\}$ is said to be
symmetric if all matrices $A_{i}$ are symmetric. Two association schemes
$\mathscr{A}$ and $\mathscr{B}$ are said to be isomorphic if there exists a
permutation matrix $P$ such that $\mathscr{B}=P^{-1}\mathscr{A}P$, otherwise,
they are called non-isomorphic.

For any $n\equiv 0\pmod{4}$, there is a known correspondence between skew
Hadamard matrices of order $n$ and nonsymmetric association schemes of order
$n-1$ and class $2$ [8]. Let $H$ be a skew Hadamard matrix $H$ of order $n$.
Denote by $D_{i}$ the diagonal matrix whose diagonal entries are the $i$th row
of $H$ and denote by $R_{i}$ the $(n-1)\times(n-1)$ matrix obtained from
$D_{i}^{-1}HD_{i}$ by deleting the $i$th row and column. Define

$A_{0}=I,\,A_{1}=\frac{1}{2}(J-2I+R_{i})\,\text{ and
}\,A_{2}=\frac{1}{2}(J-R_{i}).$

Then, $\mathscr{A}(H,i)=\\{A_{0},A_{1},A_{2}\\}$ is an association scheme of
order $n-1$ and class $2$.

Let $H_{1},\ldots,H_{157132}$ denote the $157132$ skew Hadamard matrices of
order $36$ given in Table 2. Then,

$\mathscr{A}(35)=\big{\\{}\mathscr{A}(H_{i},j)\,\big{|}\,i\in\\{1,\ldots,157132\\}\,\text{
and }\,j\in\\{1,\ldots,36\\}\big{\\}}.$

is a set of nonsymmetric association schemes of order $35$ and class $2$. We
describe how to classify nonsymmetric association schemes
$\mathscr{A}(H_{i},j)$ constructed from
$H_{i}\in\\{H_{1},\ldots,H_{157132}\\}$. Let
$\mathscr{A}=\\{A_{0},A_{1},A_{2}\\}$ and
$\mathscr{B}=\\{B_{0},B_{1},B_{2}\\}$ be nonsymmetric association schemes of
order $35$ and class $2$. Let $\mathrm{dg}(M)$ denote the digraph
corresponding to a $35\times 35$ matrix $M$ with entries in $\\{0,1\\}$. Then,
$\mathscr{A}$ and $\mathscr{B}$ are isomorphic if and only if
$\mathrm{dg}(A_{1})\simeq\mathrm{dg}(B_{1})$ or
$\mathrm{dg}(A_{1})\simeq\mathrm{dg}(B_{2})$, where $X\simeq Y$ means that two
digraphs $X$ and $Y$ are isomorphic [8]. Let $\mathrm{clg}(A)$ denote the
canonically labeled graph corresponding to $\mathrm{dg}(A)$, see [10, 17] for
the definition. It is known that $\mathrm{clg}(A)=\mathrm{clg}(B)$ if and only
if $\mathrm{dg}(A)\simeq\mathrm{dg}(B)$ [10, 17]. Then, $\mathscr{A}$ and
$\mathscr{B}$ are isomorphic if and only if
$\\{\mathrm{clg}(A_{1}),\mathrm{clg}(A_{2})\\}=\\{\mathrm{clg}(B_{1}),\mathrm{clg}(B_{2})\\}$.
Considering canonically labeled graphs significantly reduces the amount of our
computation. In this way, by comparing canonically labeled graphs
$\\{\mathrm{clg}(A_{1}),\mathrm{clg}(A_{2})\\}$ in
$\mathscr{A}=\\{A_{0},A_{1},A_{2}\\}$, we determined by Magma [2] that
$\mathscr{A}(35)$ contains $2793032$ non-isomorphic nonsymmetric association
schemes. This calculation was also done by using nauty [17].

######  Proposition 6.1.

At least $2793032$ non-isomorphic nonsymmetric association schemes of order
$35$ and class $2$ exist.

The obtained list of the non-isomorphic nonsymmetric association schemes of
order $35$ and class $2$ is available electronically at [9].

## Acknowledgments

Masaaki Harada was supported by JSPS KAKENHI under grant numbers 19H01802 and
23H01087. Hadi Kharaghani was supported by the Natural Sciences and
Engineering Research Council of Canada (NSERC). Ali Mohammadian was supported
by the Natural Science Foundation of Anhui Province with grant identifier
2008085MA03 and by the National Natural Science Foundation of China with grant
number 12171002.

## References

  * [1] M. Araya and M. Harada, Some restrictions on the weight enumerators of near-extremal ternary self-dual codes and quaternary Hermitian self-dual codes, Des. Codes Cryptogr. 91 (2023), 1813–1843.
  * [2] W. Bosma, J. Cannon and C. Playoust, The Magma algebra system I: The user language, J. Symbolic Comput. 24 (1997), 235–265.
  * [3] E. Brown and K.B. Reid, Doubly regular tournaments are equivalent to skew Hadamard matrices, J. Combin. Theory Ser. A 12 (1972), 332–338.
  * [4] J. Cooper, J. Milas and W.D. Wallis, Hadamard equivalence, Combinatorial Mathematics, in: Lecture Notes in Math., vol. 686, Springer, pp. 126–135, 1978.
  * [5] D.Ž. Ðoković, Skew-Hadamard matrices of order $276$, available at: http://arxiv.org/pdf/2301.02751.pdf.
  * [6] P. Erdős and L. Moser, On the representation of directed graphs as unions of orderings, Magyar Tud. Akad. Mat. Kutató Int. Közl. 9 (1964), 125–132.
  * [7] I.A. Faradžev, Constructive enumeration of combinatorial objects, in: Problèmes combinatoires et théorie des graphes (Colloq. Internat. CNRS, Univ. Orsay, Orsay, 1976), Colloq. Internat. CNRS, 260, CNRS, Paris, pp. 131–135, 1978.
  * [8] A. Hanaki, H. Kharaghani, A. Mohammadian and B. Tayfeh-Rezaie, Classification of skew-Hadamard matrices of order $32$ and association schemes of order $31$, J. Combin. Des. 28 (2020), 421–427.
  * [9] M. Harada, Codes and association schemes, available at: http://www.math.is.tohoku.ac.jp/~mharada/F3-H36/.
  * [10] P. Kaski and P.R.J. Östergård, Classification Algorithms for Codes and Designs, in; Algorithms Comput. Math., 15, Springer, Berlin, 2006.
  * [11] H. Kharaghani and B. Tayfeh-Rezaie, Hadamard matrices of order $32$, J. Combin. Des. 21 (2013), 212–221.
  * [12] H. Kharaghani and B. Tayfeh-Rezaie, On the classification of Hadamard matrices of order $32$, J. Combin. Des. 18 (2010), 328–336.
  * [13] H. Kharaghani and B. Tayfeh-Rezaie, A Hadamard matrix of order $428$, J. Combin. Des. 13 (2005), 435–440.
  * [14] H. Kimura, Classification of Hadamard matrices of order $28$, Discrete Math. 133 (1994), 171–180.
  * [15] E.G. Koçer, Circulant, negacyclic and semicirculant matrices with the modified Pell, Jacobsthal and Jacobsthal-Lucas numbers, Hacet. J. Math. Stat. 36 (2007), 133–142.
  * [16] C.L. Mallows and N.J.A. Sloane, An upper bound for self-dual codes, Inform. and Comput. 22 (1973), 188–200.
  * [17] B.D. McKay and A. Piperno, Practical graph isomorphism, II, J. Symbolic Comput. 60 (2014), 94–112.
  * [18] T. Miezaki, A. Munemasa and H. Nakasora, A note on Assmus–Mattson type theorems, Des. Codes Cryptogr. 89 (2021), 843–858.
  * [19] A. Mohammadian and B. Tayfeh-Rezaie, Hadamard matrices with few distinct types, Linear Multilinear Algebra 67 (2019), 1596–1605.
  * [20] R.E.A.C. Paley, On orthogonal matrices, J. Math. Phys. 12 (1933), 311–320.
  * [21] E.M. Rains and N.J.A. Sloane, Self-dual codes, in: Handbook of coding theory, Vol. I, Elsevier, Amsterdam, pp. 177–294, 1998.
  * [22] R.C. Read, Every one a winner or how to avoid isomorphism search when cataloguing combinatorial configurations, Ann. Discrete Math. 2 (1978), 107–120.
  * [23] J. Wallis, Some $(1,-1)$ matrices, J. Combin. Theory Ser. B 10 (1971), 1–11.


# Solution of nonlinear system of equations through homotopy path

A. Duttaa,1 and A. K. Dasb,2
_a Department of Mathematics, Jadavpur University, Kolkata, 700 032, India_
_b SQC & OR Unit, Indian Statistical Institute, Kolkata, 700 108, India_
_1 Email<EMAIL_ADDRESS>
_2 Email<EMAIL_ADDRESS>

###### Abstract

The paper aims to show the equivalency between nonlinear complementarity
problem and the system of nonlinear equations. We propose a homotopy method
with vector parameter $\lambda$ in finding the solution of nonlinear
complementarity problem through a system of nonlinear equations. We propose a
smooth and bounded homotopy path to obtain solution of the system of nonlinear
equations under some conditions. An oligopolistic market equilibrium problem
is considered to show the effectiveness of the proposed homotopy continuation
method. Keywords: Nonlinear complementarity problem, system of nonlinear
equations, homotopy function with vector parameter, bounded smooth curve,
oligopolistic market equilibrium.

AMS subject classifications: 90C33, 90C30, 15A69.

##  1 Introduction

A class of problems arising in various fields of sciences, can be studied via
the nonlinear system of equations using various techniques. In recent years,
researchers are interested to solve nonlinear system of equations both
analytically and numerically. Several iterative methods have been developed
using different techniques such as Taylor’s series expansion, quadrature
formulas, homotopy method, interpolation, decomposition and its various
modification. For details, see [30], [28], [8], [29], [7] and [34].

Eaves and Saigal [6] formed an important class of globally convergent methods
for solving systems of non-linear equations, which is known as homotopy
method. Such methods have been used to constructively prove the existence of
solutions to many economic and engineering problems. Let $Y,X$ be two
topologocal spaces and $p,q:Y\to X$ be continuous maps. A homotopy from $p$ to
$q$ is a continuous function $H:Y\times[0,1]$ $\to X$ satisfying
$H(y,0)=p(y),$ $H(y,1)=q(y)\ \forall\ y\in Y.$ If such a homotopy exists, then
$p$ is homotopic to $q$ and it is denoted by $p\simeq q.$ Let $p,q:R\to R$ any
two continuous, real functions, then $p\simeq q.$ Now we define a function
$H:R\times[0,1]\to R$ by $H(y,u)=(1-u)f(y)+ug(y).$ Clearly $H$ is continuous
and $H(y,0)=p(y),$ $H(y,1)=q(y).$ Thus $H$ is a homotopy between $p$ and $q.$
Let $Y,X$ be two topological spaces and Map$(Y,X)$ be the set of all
continuous maps from $Y$ to $X.$ Homotopy is an equivalence relation on
Map$(Y,X).$ The fundamental idea of the homotopy continuation method is to
solve a problem by tracing a certain continuous path that leads to a solution
to the problem. Thus, defining a homotopy mapping that yields a finite
continuation path plays an essential role in a homotopy continuation method.
The homotopy method [32] is itself an important class of globally convergent
methods. Many homotopy methods are proposed for constructive proof of the
existence of solutions to systems of nonlinear equations, nonlinear
optimization problems, Brouwer fixed point problems, nonlinear programming,
game problem and complementarity problems [33].

We are interested in solving the nonlinear system of equations. The nonlinear
complementarity problem, is identified as an important mathematical
programming problem can be converted into nonlinear system of equations.The
idea of nonlinear complementarity problem is based on the concept of linear
complementarity problem. For recent study on this problem and applications see
[4], [25], [18], [20] and references therein. For details of several matrix
classes in complementarity theory, see [9], [10], [24], [22], [13], [12],
[19], [23], [5] and references cited therein. The problem of computing the
value vector and optimal stationary strategies for structured stochastic games
for discounted and undiscounded zero-sum games and quadratic Multi-objective
programming problem are formulated as linear complementary problems. For
details see [14], [15], [21] and [17]. The complementarity problems are
considered with respect to principal pivot transforms and pivotal method to
its solution point of view. For details see [2], [26], [3] and [27].

In this paper, we consider the modified homotopy continuation method to find
the solution of the nonlinear system of equations with vector parameter. The
paper is organized as follows. Section 2 presents some basic notations and
results. In section 3, we prove the necessary and sufficient condition, which
is a system of nonlinear equations to find the solution of nonlinear
complementarity problem. A new homotopy approach with vector parameter is
proposed to find the solution of nonlinear system of equations. In the last
section, we consider an oligopoly market equilibrium problem [16] is used to
illustrate the effectiveness of the proposed homotopy approach with vector
parameter.

##  2 Preliminaries

Consider a function $f:R^{n}\rightarrow R^{n}$ , and a vector $z\in R^{n}$
such that $f=\left[\begin{array}[]{c}f_{1}\\\ f_{2}\\\ \vdots\\\ f_{n}\\\
\end{array}\right]$ and $z=\left[\begin{array}[]{c}z_{1}\\\ z_{2}\\\ \vdots\\\
z_{n}\\\ \end{array}\right].$ The complementarity problem is to find a vector
$z\in R^{n}$ such that

$z^{T}f(z)=0,\ \ \ f(z)\geq 0,\ \ \ z\geq 0.$ (2.1)

When the function $f$ is a nonlinear function, then it is called nonlinear
complementarity problem. For details see [11].

The basic idea of homotopy method is to construct a homotopy continuation path
from the auxiliary mapping g to the object mapping f. Suppose the given
problem is to find a root of the non-linear equation f(x) = 0 and suppose g(x)
= 0 is an auxiliary equation with $g(x_{0})=0$. Then the homotopy function
$H:R^{n+1}\to R^{n}$ can be defined as $H(x,\mu)=$ $(1-\mu)f(x)+\mu g(x),$
$0\leq\mu\leq 1$ is a parameter. Then we consider the homotopy equation
$H(x,\mu)=0,$ where $(x_{0},1)$ is a known solution of the homotopy equation.
Our aim is to find the solution of the equation $f(x)=0$ from the known
solution of $g(x)=0$ by solving the homotopy equation $H(x,\mu)=0$ varrying
the values of $\mu$ from $1$ to $0$.

Now we state some results.

######  Lemma 2.1:

(Generalizations of Sard’s Theorem[1]) Let $W\subset R^{n}$ be an open set and
$f:R^{n}\to R^{q}$ be smooth. We say $z\in R^{q}$ is a regular value for $f$
if $\text{Range}Df(w)=R^{q}$ $\forall\ w\in f^{-1}(z),$ where $Df(w)$ denotes
the $n\times q$ matrix of partial derivatives of $f(w).$

######  Lemma 2.2:

(Parameterized Sard Theorem [31]) Let $W\subset R^{p},Z\subset R^{q}$ be open
sets, and let $\phi:W\times Z\to R^{m}$ be a $C^{\omega}$ mapping, where
$\omega>\text{max}\\{0,q-m\\}.$ If $0\in R^{m}$ is a regular value of $\phi,$
then for almost all $a\in W,0$ is a regular value of $\phi_{w}=\phi(w,.).$

######  Lemma 2.3:

(The inverse image theorem [31]) Let $\phi:W\subset R^{p}\to R^{q}$ be
$C^{\omega}$ mapping, where $\omega>\text{max}\\{0,p-q\\}.$ Then
$\phi^{-1}(0)$ consists of some $(p-q)$-dimensional $C^{\omega}$ manifolds.

######  Lemma 2.4:

(Classification theorem of one-dimensional smooth manifold [35]) One-
dimensional smooth manifold is diffeomorphic to a unit circle or a unit
interval.

Now in the next section we will solve the nonlinear complementarity problem
2.1 using homotopy method with vector parameter.

##  3 Main results

######  Theorem 3.1:

Let $\phi:R\to R$ be any increasing function such that $\phi(0)=0$. Then $z$
solves the complementarity problem 2.1 if and only if

$\phi((f_{i}(z)-z_{i})^{2})-\phi(f_{i}(z)|f_{i}(z)|)-\phi(z_{i}|z_{i}|)=0$
(3.2)

###### Proof.

Necessery. For each $i=1,2,\ldots,n$, either $z_{i}=0,f_{i}(z)\geq 0$ or
$f_{i}(z)=0,z_{i}\geq 0$.
If $z_{i}=0,f_{i}(z)\geq 0$ then
$\phi((f_{i}(z)-z_{i})^{2})-\phi(f_{i}(z)|f_{i}(z)|)-\phi(z_{i}|z_{i}|)$
$=\phi((f_{i}(z))^{2})-\phi(f_{i}(z)|f_{i}(z)|)=\phi((f_{i}(z))^{2})-\phi((f_{i}(z))^{2})=0$.
If $f_{i}(z)=0,z_{i}\geq 0$ then
$\phi((f_{i}(z)-z_{i})^{2})-\phi(f_{i}(z)|f_{i}(z)|)-\phi(z_{i}|z_{i}|)$
$=\phi((z_{i})^{2})-\phi(z_{i}|z_{i}|)=\phi((z_{i})^{2})-\phi((z_{i})^{2})=0$.
So the solution of $(1.1)$ satisfies $(1.2)$.
Sufficient. (a) To show that $f(z)\geq 0$ assume the contrary, i.e.
$f_{i}(z)<0$ for some $i=1,2,\dots,n$. Then

$0\leq\phi((f_{i}(z)-z_{i})^{2})$
$=\phi(f_{i}(z)|f_{i}(z)|)+\phi(z_{i}|z_{i}|)$
$=\phi(-f_{i}(z)^{2})+\phi(z_{i}|z_{i}|)$ $<\phi(z_{i}|z_{i}|)$

This implies that $\phi(z_{i}|z_{i}|)>0$ $\implies z_{i}|z_{i}|>0$ $\implies
z_{i}>0$ and $\phi((f_{i}(z)-z_{i})^{2})<\phi(z_{i}|z_{i}|)$
$\implies((f_{i}(z)-z_{i})^{2})<z_{i}|z_{i}|$.
But for $z_{i}>0$, $(z_{i})^{2}<(f_{i}(z)-z_{i})^{2}$ [as
$(f_{i}(z)-z_{i})<0,z_{i}>0$ then $(f_{i}(z)-z_{i})<0$,
$|f_{i}(z)-z_{i}|>z_{i}$ so $(f_{i}(z)-z_{i})^{2}>z_{i}^{2}$ ].
So, it contradicts that $f_{i}(z)<0$. So $f_{i}(z)\geq 0$.
(b) To show that $z\geq 0$ interchange the roles of $z_{i}$ and $f_{i}(z)$.
(c)From (a) and (b) we have that $z\geq 0$ and $f(z)\geq 0$. To show
$z^{T}f(z)=0$ assume the contrary $z_{i}>0$ and $f_{i}(z)>0$ for some
$i=1,2,\ldots,n.$
If $f_{i}(z)\geq z_{i}$, then
$\phi((f_{i}(z)-z_{i})^{2})<\phi((f_{i}(z))^{2})+\phi((z_{i})^{2})=\phi(f_{i}(z)|f_{i}(z)|)+\phi(z_{i}|z_{i}|)$.
This contradicts that
$\phi((f_{i}(z)-z_{i})^{2})=\phi(f_{i}(z)|f_{i}(z)|)+\phi(z_{i}|z_{i}|)$.
Similarly, for $z_{i}(z)\geq f_{i}(z),$ interchanging the roles of $z_{i}$ and
$f_{i}(z)$, we get a contradiction.

Hence it is shown that the complementarity problem of finding a $z\in R^{n}$
satisfying $z^{T}f(z)=0,\ \ \ f(z)\geq 0,\ \ \ z\geq 0$ where
$f:R^{n}\rightarrow R^{n}$ is equivalent to the problem of solving system of n
nonlinear equations in n variables.

${}\psi_{i}(z)=\phi((f_{i}(z)-z_{i})^{2})-\phi(f_{i}(z)|f_{i}(z)|)-\phi(z_{i}|z_{i}|)=0\
\ \ \forall i\in\\{1,2,\cdots n\\},$ (3.3)

where $\phi$ is an increasing function defined by $\phi:R\to R$ such that
$\phi(x)=x^{3}.$ Now $\frac{\partial\psi_{i}}{\partial z_{j}}$
$=\phi^{{}^{\prime}}((f_{i}(z)-z_{i})^{2})2(f_{i}(z)-z_{i})(\frac{\partial
f_{i}}{\partial z_{j}}-\delta_{ij})$
$-\phi^{{}^{\prime}}(f_{i}(z)|f_{i}(z)|)2f_{i}(z)sgn(f_{i}(z))\frac{\partial
f_{i}}{\partial z_{j}}$
$-\phi^{{}^{\prime}}(z_{i}|z_{i}|)2z_{i}sgn(z_{i})\delta_{ij}$ where

$sgn(x)=\begin{cases}\ 1&\quad\text{if }x>0\\\ \ -1&\quad\text{if }x<0\\\
\end{cases}$

######  Theorem 3.2:

Let $z$ be the nondegenerate solution of nonlinear complementarity problem
2.1, i.e. $z+f(z)>0$. Let $\cal{J}$ $(f(z))$, the Jacobian of $f$ at $z$, have
nonsingular principal minors and let $\phi:R\to R$ be a differentiable
strictly increasing function such that $\phi^{\prime}(y)>0$ for all $y>0$ and
$\phi(0)=0.$ Then $z$ solves 3.3 and the Jacobian of $\psi$ at $z$, $\cal{J}$
$(\psi(z))$ is nonsingular.

###### Proof.

The Jacobian matrix of $\psi(z)$ is defined by $\cal{J}$ $(\psi(z))=$
$\left[\begin{array}[]{rrrr}\frac{\partial\psi_{1}}{\partial
z_{1}}&\frac{\partial\psi_{1}}{\partial
z_{2}}&\cdots&\frac{\partial\psi_{1}}{\partial z_{n}}\\\
\frac{\partial\psi_{2}}{\partial z_{1}}&\frac{\partial\psi_{2}}{\partial
z_{2}}&\cdots&\frac{\partial\psi_{2}}{\partial z_{n}}\\\
\vdots&\vdots&\vdots&\vdots\\\ \frac{\partial\psi_{n}}{\partial
z_{1}}&\frac{\partial\psi_{n}}{\partial
z_{2}}&\cdots&\frac{\partial\psi_{n}}{\partial z_{n}}\end{array}\right],$

where $\frac{\partial\psi_{i}}{\partial
z_{i}}=\phi^{{}^{\prime}}((f_{i}(z)-z_{i})^{2})2(f_{i}(z)-z_{i})(\frac{\partial
f_{i}}{\partial
z_{i}}-1)-\phi^{{}^{\prime}}(f_{i}(z)|f_{i}(z)|)2f_{i}(z)sgn(f_{i}(z))\frac{\partial
f_{i}}{\partial z_{i}}-\phi^{{}^{\prime}}(z_{i}|z_{i}|)2z_{i}sgn(z_{i})$ and
$\frac{\partial\psi_{i}}{\partial
z_{j}}=\phi^{{}^{\prime}}((f_{i}(z)-z_{i})^{2})2(f_{i}(z)-z_{i})(\frac{\partial
f_{i}}{\partial z_{j}})-$
$\phi^{{}^{\prime}}(f_{i}(z)|f_{i}(z)|)2f_{i}(z)sgn(f_{i}(z))\frac{\partial
f_{i}}{\partial z_{j}}$ for $i\neq j.$
Assume that $f_{i}(z)=0$ for $i=1,2,\cdots,n_{1},n_{1}\leq n$ and $f_{i}(z)>0$
for $i=n_{1}+1,n_{1}+2,\cdots n.$ Hence by nondegeneracy of $z$, $z_{i}>0$ for
$i=1,2,\cdots n_{1},n_{1}\leq n,$ $z_{i}=0$ for $i=n_{1}+1,n_{1}+2,\cdots n.$
Now $\cal{J}$ $(\psi(z))=$
$\left[\begin{array}[]{rrrr}-\phi^{{}^{\prime}}((z_{1})^{2})2z_{1}&0&\cdots&0\\\
0&-\phi^{{}^{\prime}}((z_{2})^{2})2z_{2}&\cdots&0\\\
\vdots&\vdots&\cdots&\vdots\\\
0&0&\cdots&-\phi^{{}^{\prime}}((z_{n_{1}})^{2})2z_{n_{1}}\\\ 0&0&\cdots&0\\\
\vdots&\vdots&\cdots&\vdots\\\ 0&0&\cdots&0\\\ \end{array}\right]$ $\cal{J}$
$(f(z))$ $+$ $\left[\begin{array}[]{rrrr}0&0&\cdots&0\\\
\vdots&\vdots&\cdots&\vdots\\\ 0&0&\cdots&0\\\
-\phi^{{}^{\prime}}(({f_{n_{1}+1}}(z))^{2})2{f_{n_{1}+1}}(z)&0&\cdots&0\\\
0&-\phi^{{}^{\prime}}(({f_{n_{1}+2}}(z))^{2})2{f_{n_{1}+2}}(z)&\cdots&0\\\
\vdots&\vdots&\cdots&\vdots\\\
0&0&\cdots&-\phi^{{}^{\prime}}(({f_{n}}(z))^{2})2{f_{n}}(z)\\\ 0&0&\cdots&0\\\
\vdots&\vdots&\cdots&\vdots\\\ 0&0&\cdots&0\\\ \end{array}\right]$. Since
$\cal{J}$ $(f(z))$ has nonsingular principal minors and $\phi^{\prime}(y)>0$
for $y>0,$ the Jacobian of $\psi(z),$ $\cal{J}$ $(\psi(z))$ is nonsingular.

###  3.1 Homotopy based solution with vector parameter $\tilde{\lambda}\in
R^{n}$

The basic idea of homotopy method with $\tilde{\lambda}\in R^{n}$ is to
construct a multidimentional homotopy path to find the solution of the object
fuction $\tilde{f}(\tilde{x})=0$ varying each component of $\tilde{\lambda}$
from $1$ to $0.$ We consider the homotopy function
$\mathcal{H}(\tilde{x},\tilde{x}^{(0)},\tilde{\lambda})=\tilde{f}(\tilde{x})-\tilde{\lambda}\tilde{f}(\tilde{x}^{(0)}),$
where each component of $\tilde{f}(\tilde{x}^{(0)}),\ \
f_{i}(\tilde{x}^{(0)})\neq 0\ \forall i\in\\{1,2,\cdots,n\\}$ and the product
term $\tilde{\lambda}\tilde{f}(\tilde{x}^{(0)})$ is a componentwise product
that is
$\tilde{\lambda}\tilde{f}(\tilde{x}^{(0)})=\left[\begin{array}[]{c}{\lambda_{1}}f_{1}(\tilde{x}^{(0)})\\\
{\lambda_{2}}f_{2}(\tilde{x}^{(0)})\\\ \vdots\\\
{\lambda_{n}}f_{n}(\tilde{x}^{(0)})\\\ \end{array}\right],$ where
$\tilde{\lambda}=\left[\begin{array}[]{c}{\lambda_{1}}\\\ {\lambda_{2}}\\\
\vdots\\\ {\lambda_{n}}\\\ \end{array}\right]$ and
$\tilde{f}(\tilde{x}^{(0)})=\left[\begin{array}[]{c}f_{1}(\tilde{x}^{(0)})\\\
f_{2}(\tilde{x}^{(0)})\\\ \vdots\\\ f_{n}(\tilde{x}^{(0)})\\\
\end{array}\right].$
In this method our main aim is to vary each component
$\lambda_{i}=\frac{f_{i}(\tilde{x})}{f_{i}(\tilde{x}^{(0)})}$ of
$\tilde{\lambda}$ from $1$ to $0.$

Now we solve the nonlinear equation 3.2 using the homotopy function with
vector parameter
$\tilde{\lambda}=[\lambda_{1},\lambda_{2},\cdots,\lambda_{n}]^{t}\in R^{n}.$
Let
$\psi_{i}(z)=\phi((f_{i}(z)-z_{i})^{2})-\phi(f_{i}(z)|f_{i}(z)|)-\phi(z_{i}|z_{i}|)=0.$
Consider the set $F_{\mathcal{H}}=\\{\tilde{x}:\tilde{x}\in R^{n}\\}$ and
$\tilde{F}_{\mathcal{H}}=\\{(\tilde{x},\tilde{\lambda}):(\tilde{x},\tilde{\lambda})\in
R^{n}\times(0,1]^{n}\\}.$ The homotopy function can be given as,

${}\mathcal{H}(\tilde{x},\tilde{x}^{(0)},\tilde{\lambda})=\psi(\tilde{x})-\tilde{\lambda}\psi(\tilde{x}^{(0)})=\left[\begin{array}[]{c}\psi_{1}(\tilde{x})-{\lambda_{1}}\psi_{1}(\tilde{x}^{(0)})\\\
\psi_{2}(\tilde{x})-{\lambda_{2}}\psi_{2}(\tilde{x}^{(0)})\\\ \vdots\\\
\psi_{n}(\tilde{x})-{\lambda_{n}}\psi_{n}(\tilde{x}^{(0)})\end{array}\right]$
(3.4)

Where, $\tilde{x}^{(0)}$ is the initial known value of the vector
$\tilde{x}\in R^{n}$ such that $\psi_{i}(\tilde{x}^{(0)})\neq 0\ \forall i\
\in\\{1,2,\cdots,n\\}$ and $\tilde{\lambda}\in R^{n}.$ In this modification,
$\tilde{\lambda}$ is taken as a vector.

######  Theorem 3.3:

If the Jacobian matrix
$\frac{\partial{\psi}(\tilde{x}^{(0)})}{\partial\tilde{x}^{(0)}}$ at initial
point $\tilde{x}^{(0)}$ is nonsingular, then for almost all initial points
$\tilde{x}^{(0)}\in\tilde{F}_{\mathcal{H}},$ $0$ is a regular value of the
homotopy function $\mathcal{H}:R^{6n}\times(0,1]^{n}\to R^{6n}$ and the zero
point set
$\mathcal{H}_{\tilde{x}^{(0)}}^{-1}(0)=\\{(\tilde{x},\tilde{\lambda})\in\tilde{F}_{\mathcal{H}}:\mathcal{H}_{\tilde{x}^{(0)}}(\tilde{x},\tilde{\lambda})=0\\}$
contains a smooth curve $\Gamma_{\tilde{x}^{(0)}}$ starting from
$(\tilde{x}^{(0)},e).$

###### Proof.

The Jacobian matrix of the above homotopy function3.4
$\mathcal{H}(\tilde{x},\tilde{x}^{(0)},\tilde{\lambda})$ is denoted by
$D\mathcal{H}(\tilde{x},\tilde{x}^{(0)},\tilde{\lambda}))$ and we have
$D\mathcal{H}(\tilde{x},\tilde{x}^{(0)},\tilde{\lambda}))=$
$\left[\begin{array}[]{ccc}\frac{\partial{\mathcal{H}(\tilde{x},\tilde{x}^{(0)},\tilde{\lambda})}}{\partial{\tilde{x}}}&\frac{\partial{\mathcal{H}(\tilde{x},\tilde{x}^{(0)},\tilde{\lambda})}}{\partial{\tilde{x}^{(0)}}}&\frac{\partial{\mathcal{H}(\tilde{x},\tilde{x}^{(0)},\tilde{\lambda})}}{\partial{\tilde{\lambda}}}\\\
\end{array}\right].$
For all $\tilde{x}^{(0)}\in\tilde{F}_{H}$ such that $\psi(\tilde{x}^{(0)})\neq
0$ and $\tilde{\lambda}\in(0,1]^{n},$ we have
$\frac{\partial{\mathcal{H}(\tilde{x},\tilde{x}^{(0)},\tilde{\lambda})}}{\partial{\tilde{x}^{(0)}}}=$
$\left[\begin{array}[]{cccc}-{\lambda_{1}}\frac{\partial{\psi_{1}}(\tilde{x}^{(0)})}{\partial\tilde{x}_{1}^{(0)}}\
-{\lambda_{1}}\frac{\partial{\psi_{1}}(\tilde{x}^{(0)})}{\partial\tilde{x}_{2}^{(0)}}\
\cdots\
-{\lambda_{1}}\frac{\partial{\psi_{1}}(\tilde{x}^{(0)})}{\partial\tilde{x}_{n}^{(0)}}\\\
-{\lambda_{2}}\frac{\partial{\psi_{2}}(\tilde{x}^{(0)})}{\partial\tilde{x}_{1}^{(0)}}\
-{\lambda_{2}}\frac{\partial{\psi_{2}}(\tilde{x}^{(0)})}{\partial\tilde{x}_{2}^{(0)}}\
\cdots\
-{\lambda_{2}}\frac{\partial{\psi_{2}}(\tilde{x}^{(0)})}{\partial\tilde{x}_{n}^{(0)}}\\\
\vdots\ \vdots\ \vdots\ \vdots\\\
-{\lambda_{n}}\frac{\partial{\psi_{n}}(\tilde{x}^{(0)})}{\partial\tilde{x}_{1}^{(0)}}\
-{\lambda_{n}}\frac{\partial{\psi_{n}}(\tilde{x}^{(0)})}{\partial\tilde{x}_{2}^{(0)}}\
\cdots\
-{\lambda_{n}}\frac{\partial{\psi_{n}}(\tilde{x}^{(0)})}{\partial\tilde{x}_{n}^{(0)}}\end{array}\right].$
Now
$\det\frac{\partial{\mathcal{H}(\tilde{x},\tilde{x}^{(0)},\tilde{\lambda})}}{\partial{\tilde{x}^{(0)}}}=$
$(-1)^{n}\det J_{\mathcal{H}}{\tilde{x}^{(0)}}\prod_{i=1}^{i=n}\lambda_{i}\neq
0$ for $\lambda\in(0;1]^{n},$ where
$J_{\mathcal{H}}{\tilde{x}^{(0)}}=\frac{\partial{\psi}(\tilde{x}^{(0)})}{\partial\tilde{x}^{(0)}}.$
Thus $D\mathcal{H}(\tilde{x},\tilde{x}^{(0)},\tilde{\lambda}))$ is of full row
rank. Therefore, $0$ is a regular value of
$\mathcal{H}(\tilde{x},\tilde{x}^{(0)},\tilde{\lambda})$ and
$\mathcal{H}_{\tilde{x}^{(0)}}^{-1}(0)$ consists of some smooth curves such
that $\mathcal{H}_{\tilde{x}^{(0)}}(\tilde{x}^{(0)},1)=0.$ Hence there must be
a smooth curve $\Gamma_{\tilde{x}^{(0)}}$ starting from $(\tilde{x}^{(0)},1).$

######  Theorem 3.4:

Assume that the function $f(\tilde{x})$ is either increasing or bounded
function. The smooth curve $\Gamma_{\tilde{x}^{(0)}}$ is a bounded curve.

###### Proof.

From the homotopy function 3.4 we get
$\psi_{i}(\tilde{x})={\lambda_{i}}\psi_{i}(\tilde{x}^{(0)})\ \forall
i\in\\{1,2,\cdots,n\\},$where
$\psi_{i}(\tilde{x})=\phi((f_{i}(\tilde{x})-\tilde{x}_{i})^{2})-\phi(f_{i}(\tilde{x})|f_{i}(\tilde{x})|)-\phi(\tilde{x}_{i}|\tilde{x}_{i}|).$
It is clear that $\lambda_{i}\in[0,1],$ $\psi_{i}(\tilde{x})$ is finite. So
$\|\psi(\tilde{x})\|<\infty.$
This implies that
$\psi_{i}(\tilde{x})=\phi((f_{i}(\tilde{x})-\tilde{x}_{i})^{2})-\phi(f_{i}(\tilde{x})|f_{i}(\tilde{x})|)-\phi(\tilde{x}_{i}|\tilde{x}_{i}|)={\lambda_{i}}\psi_{i}(\tilde{x}^{(0)})\
\forall i\in\\{1,2,\cdots,n\\}.$ Let
${\lambda_{i}}\psi_{i}(\tilde{x}^{(0)})=k_{i}.$ For
$\lambda_{i}\in(0,1],k_{i}>0.$
Hence
$\phi((f_{i}(\tilde{x})-\tilde{x}_{i})^{2})=\phi(f_{i}(\tilde{x})|f_{i}(\tilde{x})|)+\phi(\tilde{x}_{i}|\tilde{x}_{i}|)+k_{i}.$
Again $\forall i\in\\{1,2,\cdots,n\\}$ we have
$\psi_{i}(\tilde{x})=((f_{i}(\tilde{x})-\tilde{x}_{i})^{6})-(f_{i}(\tilde{x})|f_{i}(\tilde{x})|)^{3}-(\tilde{x}_{i}|\tilde{x}_{i}|)^{3}=(f_{i}(\tilde{x}))^{6}-6(f_{i}(\tilde{x}))^{5}\tilde{x}_{i}+15(f_{i}(\tilde{x}))^{4}(\tilde{x}_{i})^{2}-20(f_{i}(\tilde{x}))^{3}(\tilde{x}_{i})^{3}+15(f_{i}(\tilde{x}))^{2}(\tilde{x}_{i})^{4}-6(f_{i}(\tilde{x}))(\tilde{x}_{i})^{5}+(\tilde{x}_{i})^{6}-(f_{i}(\tilde{x}))^{5}|f_{i}(\tilde{x})|-(\tilde{x}_{i})^{5}|(\tilde{x}_{i})|.$
Consider the homotopy path $\Gamma_{\tilde{x}^{(0)}}$ is unbounded. Then there
exists a sequence of points
$(\tilde{x}^{(k)},\tilde{\lambda}^{(k)})\subset\Gamma_{\tilde{x}^{(0)}}$ such
that $\|\tilde{x}^{(k)}\|\to\infty.$ Then there are two possibilities.
Case 1: Let $\|\tilde{x}^{(k)}\|\to\infty.$ Then $\exists\
i\in\\{1,2,\cdots,n\\}$ such that $\tilde{x}^{(k)}_{i}\to-\infty$ as
$k\to\infty.$ Let set
$I_{1\tilde{x}}=\\{i\in\\{1,2,\cdots,n\\}:\lim\limits_{k\to\infty}\tilde{x}^{(k)}_{i}\to-\infty\\}.$
Now consider $\|f(\tilde{x}^{(k)})\|<\infty.$ Then from above we have $\forall
i\in I_{1\tilde{x}},$ as $k\to\infty,$
$\frac{\psi_{i}(\tilde{x}^{(k)})}{(\tilde{x}^{(k)}_{i})^{6}}\to 2,$ which
contradicts that $\psi(x)$ is bounded.
Again consider that the nonlinear function $f(\tilde{x})$ is unbounded. It is
noted that $f(\tilde{x})$ is an increasing function.
Let
$L_{1f(\tilde{x})}=\\{l\in\\{1,2,\cdots,n\\}:\lim\limits_{k\to\infty}f_{l}(\tilde{x}^{(k)})\to-\infty\\}$
and consider
$\lim\limits_{k\to\infty}\frac{f_{i}(\tilde{x}^{(k)})}{(\tilde{x}_{i}^{(k)})}=p\
\forall\ i\in I_{1\tilde{x}}\cap L_{1f(\tilde{x})}$. Therefore $\forall\ i\in
I_{1\tilde{x}}\cap L_{1f(\tilde{x})},$ as $k\to\infty,$
$\frac{\psi_{i}(\tilde{x}^{(k)})}{(\tilde{x}^{(k)}_{i})^{6}}\to
2p^{6}-6p^{5}+15p^{4}-20p^{3}+15p^{2}-6p+2.$ From the boundedness of
$\psi(x^{(k)})$, it is clear that
$2p^{6}-6p^{5}+15p^{4}-20p^{3}+15p^{2}-6p+2=0,$ has no real solution, which is
a contradiction.
Now for all $l\in L_{1f(\tilde{x})},l\notin I_{1\tilde{x}},$
$\lim\limits_{k\to\infty}\frac{\psi_{l}(\tilde{x}^{(k)})}{f_{l}(\tilde{x}^{(k)})^{6}}\to
2,$ contradicts the boundedness of the function $\psi(\tilde{x}).$

Case 2: Let $\|\tilde{x}^{(k)}\|\to\infty.$ Then $\exists\
j\in\\{1,2,\cdots,n\\}$ such that $\tilde{x}^{(k)}_{j}\to\infty$ as
$k\to\infty.$ Let set
$I_{2\tilde{x}}=\\{j\in\\{1,2,\cdots,n\\}:\lim\limits_{k\to\infty}\tilde{x}^{(k)}_{j}\to\infty\\}$.
Let $\tilde{x}_{j}>0,{f_{j}}(\tilde{x})\geq 0.$ Consider that
$\lambda_{j}\in(0,1].$Then
$\phi((f_{j}(\tilde{x})-\tilde{x}_{j})^{2})=\phi(f_{j}(\tilde{x})|f_{j}(\tilde{x})|)+\phi(\tilde{x}_{j}|\tilde{x}_{j}|)+k_{j}=\phi((f_{j}(\tilde{x})^{2})+\phi((\tilde{x}_{j})^{2})+k_{j},$
where $k_{j}>0.$ Now
$0\leq\phi((f_{j}(\tilde{x})-\tilde{x}_{j})^{2})\leq\phi((f_{j}(\tilde{x})^{2})+\phi((\tilde{x}_{j})^{2}).$
This contradicts that
$\phi((f_{j}(\tilde{x})-\tilde{x}_{j})^{2})=\phi((f_{j}(\tilde{x})^{2})+\phi((\tilde{x}_{j})^{2})+k_{j}.$
Therefore if $\psi_{j}(\tilde{x})={\lambda_{j}}\psi_{j}(\tilde{x}^{(0)}),$
$\lambda_{j}\in(0,1]$ has the solution $\tilde{x}_{j}>0,$ then
$f_{j}(\tilde{x})<0.$
Now consider that $\|f(\tilde{x}^{(k)})\|<\infty.$
Then $\forall j\in I_{2\tilde{x}},$ as $k\to\infty,$
$\frac{\psi_{j}(\tilde{x}^{(k)})}{(\tilde{x}^{(k)}_{j})^{5}}\to-6(f_{j}(\tilde{x}^{(k)}))\nrightarrow
0,$ which is a contradiction.
Therefore considering both the cases we conclude that the homotopy path
$\Gamma_{\tilde{x}^{(0)}}$ is bounded for $\lambda_{i}\in(0,1]\ \forall i$.

Note that here we use the vector $\tilde{\lambda},$ such that
$\tilde{\lambda}:[0,\infty]^{n}\to(0,1]^{n}$ with $\lambda_{i}=\exp(-t_{i}).$

######  Remark 3.1:

Now we trace the homotopy path
$\Gamma_{\tilde{x}^{(0)}}\subset\mathcal{H}_{\tilde{x}^{(0)}}^{-1}(0)\subset\tilde{F}_{\mathcal{H}}$
from the initial point $(\tilde{x}^{(0)},e)$ until
$\tilde{\lambda}\to\tilde{0}$ and find the solution of the system of nonlinear
equation 3.2 under some assumptions. If the homotopy path is bounded and
$\tilde{\lambda}$ goes to $0$ starting from $e,$ the component $\tilde{x}$
give the solution of 3.2. Let $s$ denote the arc length of
$\Gamma_{\tilde{x}^{(0)}}$, we can parameterize the homotopy path
$\Gamma_{\tilde{x}^{(0)}}$ with respect to $s$ in the following form

$\mathcal{H}_{\tilde{x}^{(0)}}(\tilde{x}(s),\tilde{\lambda}(s))=0,\
\tilde{x}(0)=\tilde{x}^{(0)},\ \tilde{\lambda}(0)=e.$ (3.5)

Now differentiating 3.5 with respect to $s$, we obtain the following system of
ordinary differential equations with given initial values

$\mathcal{H}^{\prime}_{\tilde{x}^{(0)}}(\tilde{x}(s),\tilde{\lambda}(s))\left[\begin{array}[]{c}\frac{d\tilde{x}}{ds}\\\
\frac{d\tilde{\lambda}}{ds}\\\ \end{array}\right]=0,\
\|(\frac{d\tilde{x}}{ds},\frac{d\tilde{\lambda}}{ds})\|_{2}=1,\
\tilde{x}(0)=\tilde{x}^{(0)},\ \tilde{\lambda}(0)=e.\ $ (3.6)

and the $\tilde{x}$-component of
$(\tilde{x}(\bar{s}),\tilde{\lambda}(\bar{s}))$ gives the solution of the
complementarity problem for $\tilde{\lambda}({\bar{s}})=0.$

######  Theorem 3.5:

The homotopy path $\Gamma_{\tilde{x}^{(0)}}$ is determined by the initial
value problem 3.6.

###### Proof.

Let $s$ denote the arc length of the homotopy path $\Gamma_{\tilde{x}^{(0)}}.$
Now differentiating the homotopy equation 3.5, we get
$\frac{\partial\mathcal{H}}{\partial\tilde{x}}\frac{d\tilde{x}}{ds}+\frac{\partial\mathcal{H}}{\partial\tilde{\lambda}}\frac{d\tilde{\lambda}}{ds}=0.$
Let $\nu=\left[\begin{array}[]{c}\tilde{x}\\\ \tilde{\lambda}\\\
\end{array}\right].$ As $s$ is the arc length of $\Gamma_{\tilde{x}^{(0)}},$
then
$\|\nu^{\prime}\|=\|(\frac{d\tilde{x}}{ds},\frac{d\tilde{\lambda}}{ds})\|_{2}=1.$
Then we get the following system of ordinary differential equation.

$\left[\begin{array}[]{cc}\frac{\partial\mathcal{H}}{\partial\tilde{x}}&\frac{\partial\mathcal{H}}{\partial\tilde{\lambda}}\\\
\end{array}\right]\mu=0,\mu^{t}\mu=1,\frac{d\nu}{ds}=\mu,\nu(0)=\left[\begin{array}[]{c}\tilde{x}^{(0)}\\\
e\\\ \end{array}\right]$ (3.7)

Hence from the system 3.7 the first two equations are solvable on $\mu.$
Solving the following cauchy problem, the homotopy curve
$\Gamma_{\tilde{x}^{(0)}}$ can be derived.

$\frac{d\nu}{ds}=\mu,\nu(0)=\left[\begin{array}[]{c}\tilde{x}^{(0)}\\\ e\\\
\end{array}\right].$ (3.8)

###  3.2 Algorithm

Algorithm 1 Homotopy Continuation Method

Step 0: Set $i=i_{c}=0.$ [$i$ is the Number of Iteration(s) and $i_{c}$ is the
Number of shifting of the Initial Point(s).] Give an initial point
$(\tilde{x}^{(0)},\tilde{\lambda}_{0})\in\tilde{F}_{\mathcal{H}}\times\\{1\\}^{n}.$
Set $\eta_{1}=10^{-12},c_{0}=50.$

$\kappa_{1}=\sqrt{2},\kappa_{2}=9000,$ where the step-length is determined by
$\kappa_{1}^{k},k\in Z$ and the limit of the maximum step-length is maintained
by $\kappa_{1}^{k}\leq\kappa_{2}.$

Set $\tilde{u}_{1}=e,$ $\tilde{p}_{1}=\tilde{\lambda}^{-1}(e)$ and
$\tilde{p}_{0}=\tilde{\lambda}^{-1}(\tilde{0}),$ where
$e=\left[\begin{array}[]{c}1\\\ 1\\\ \vdots\\\ 1\\\ \end{array}\right],$
vector of all $1^{\prime}$s and $\tilde{0}=\left[\begin{array}[]{c}0\\\ 0\\\
\vdots\\\ 0\\\ \end{array}\right],$ vector of all $0^{\prime}$s.

Set $\epsilon_{1}=10^{-16}e,\epsilon_{2}=10^{-10}e.$ These are real numbers,
used as thresholds for $\tilde{\lambda}.$ If $\tilde{\lambda}$ achieves a
value $0\leq\tilde{\lambda}\leq\epsilon_{1},$ then the algorithm stops with an
acceptable solution. If $\tilde{\lambda}$ achieves a value, such that,
$\epsilon_{1}<\tilde{\lambda}\leq\epsilon_{2},$ the algorithm stops, declaring
that point as probable solution.

Step 1: Set $\left[\begin{array}[]{c}\tilde{x}\\\ \tilde{t}\\\
\end{array}\right]=\left[\begin{array}[]{c}\tilde{x}^{(0)}\\\ \tilde{p_{1}}\\\
\end{array}\right].$ Now calculate the constant
$d^{(0)}=\det(\frac{\partial{\mathcal{H}}}{\partial{\tilde{x}}}(\tilde{x}^{(0)},\tilde{\lambda}_{0}))$
and $s^{(0)}=\frac{d{\tilde{\lambda}}}{dt}(t=p_{1}).$ If
$|d^{(0)}|\leq\epsilon$ or $|s^{(0)}|\leq\epsilon,$ then stop else go to step
$2$ [$\epsilon\to 0,$ a threshold].

Step 2: Set $c_{1}=c_{2}=0.$ Calculate the constant
$d=\det(\frac{\partial{\mathcal{H}}}{\partial{\tilde{x}}}(\tilde{x},\tilde{\lambda})),N=\frac{d{\tilde{\lambda}(t)}}{dt},$
and $s=Ne.$ If $|d|\leq\epsilon,$ then stop else go to step $3$.

Step 3: Determine the unit predictor direction $\tau^{(n)}$ by the following
method: If sign$(d)=-\text{sign}(d_{0}),$ then set
$\tilde{t}_{d}=N^{-1}(e-\tilde{\lambda}(t)),$ else set
$\tilde{t}_{d}=-N^{-1}\tilde{\lambda}(t).$ Calculate
$\tilde{w}_{d}=-(\frac{\partial{\mathcal{H}}}{\partial{\tilde{x}}}(\tilde{x},\tilde{\lambda}))^{-1}(\frac{\partial{\mathcal{H}}}{\partial{\tilde{\lambda}}}(\tilde{x},\tilde{\lambda}))\tilde{t}_{d},$
$\tau^{(n)}=$ $\left[\begin{array}[]{c}\tilde{x}_{n}\\\ \tilde{t}_{n}\\\
\end{array}\right]=\frac{1}{\|\tilde{x}_{d},\tilde{t}_{d}\|}\left[\begin{array}[]{c}\tilde{x}_{d}\\\
\tilde{t}_{d}\\\ \end{array}\right],$
$\tilde{\tau}=\dfrac{\|\tilde{t}_{d}\|}{\|\tilde{x}_{d},\tilde{t}_{d}\|},$
where
$\|\tilde{x}_{d},\tilde{t}_{d}\|=\sqrt{\tilde{x}_{d}^{2}+\tilde{t}_{d}^{2}}.$
If $\tilde{\tau}\leq\eta_{1},$ then set $c_{1}=c_{1}+1$ else reset $c_{1}=0.$
If $c_{1}<c_{0},$ then go to step $4$ else, if $\tilde{u}_{1}\leq\epsilon_{2}$
then, stop with a probable solution else, stop due to non-convergence.

Step 4: Choosing step length: Set
$k=0;\gamma=[\nabla\mu(\tilde{x})]^{t}\tilde{x}_{n},$ where $\mu:R^{n}\to R,$
is used to increase step length in the descent direction(s). Set this
Function, such that, $\mu(\bar{\tilde{x}})\leq\mu(\tilde{x}),$ and
$\tilde{x},\bar{\tilde{x}}\in\mathcal{F}_{\mathcal{H}},$ where
$\bar{\tilde{x}}$ is the solution of the problem. $\mu(\tilde{x})$ is taken as
$[\mathcal{H}_{0}(\tilde{x})]^{t}[\mathcal{H}_{0}(\tilde{x})],$ where
$\mathcal{H}_{0}(\tilde{x})=\left[\begin{array}[]{c}\psi_{1}(\tilde{x})\\\
\psi_{2}(\tilde{x})\\\ \vdots\\\ \psi_{n}(\tilde{x})\end{array}\right].$

If $\gamma\geq 0,$
$\tilde{x}+\kappa_{1}^{k+1}\tilde{x}_{n}\in\mathcal{F}_{\mathcal{H}},\min(\tilde{p}_{0},\tilde{p}_{1})<\tilde{t}+\kappa_{1}^{k+1}\tilde{t}_{n}<\max(\tilde{p}_{0},\tilde{p}_{1}),$
then set $k=k+1$ and go to step $5.$

else if
$\gamma<0,\mu(\tilde{x}+\kappa_{1}^{k+1}\tilde{x}_{n})<\mu(\tilde{x}+\kappa_{1}^{k}\tilde{x}_{n}),\tilde{x}+\kappa_{1}^{k+1}\tilde{x}_{n}\in\mathcal{F}_{\mathcal{H}},\min(\tilde{p}_{0},\tilde{p}_{1})<\tilde{t}+\kappa_{1}^{k+1}\tilde{t}_{n}<\max(\tilde{p}_{0},\tilde{p}_{1}),$
then set $k=k+1,$ and go to step $5.$

else reset $c_{2}=0,$ and jump to step $6.$

Step 5: If $\kappa_{1}^{k}>\kappa_{2},$ then set $k=k-1,$ $c_{2}=c_{2}+1$ and
go to step $6$, else go to step $4$.

Step 6: If $c_{2}<c_{0},$ then go to step $7$, else if
$\tilde{u}_{1}\leq\epsilon_{2},$ then stop with probable solution else stop.

Step 7: Compute the predictor and corrector point:
$\left[\begin{array}[]{c}\tilde{x}_{p}\\\ \tilde{t}_{p}\\\
\end{array}\right]=\left[\begin{array}[]{c}\tilde{x}\\\ \tilde{t}\\\
\end{array}\right]+\kappa_{1}^{k}\left[\begin{array}[]{c}\tilde{x}_{n}\\\
\tilde{t}_{n}\\\ \end{array}\right],$
$\left[\begin{array}[]{c}\tilde{x}_{c}\\\ \tilde{t}_{c}\\\
\end{array}\right]=\left[\begin{array}[]{c}\tilde{x}_{p}\\\ \tilde{t}_{p}\\\
\end{array}\right]-[J_{\mathcal{H}}(\tilde{x}_{p},\tilde{t}_{p})^{+}H(\tilde{x}_{p},\tilde{t}_{p})],$
where $[J_{\mathcal{H}}(\tilde{x}_{p},\tilde{t}_{p})^{+}]$ is the Moore-
Penrose Inverse. If
$\det(\text{diag}(\mathcal{H}(\tilde{x}_{c},\tilde{p}_{0})-\mathcal{H}(\tilde{x}_{c},\tilde{p}_{1})))\leq\epsilon,$
where $\epsilon\to 0,\text{a threshold},$ then stop else
$\tilde{u}_{s}=(\text{diag}(\mathcal{H}(\tilde{x}_{c},\tilde{p}_{0})-\mathcal{H}(\tilde{x}_{c},\tilde{p}_{1}))^{-1}\mathcal{H}(\tilde{x}_{c},\tilde{p}_{0}).$
Now if $(\tilde{x}_{c},\tilde{u}_{s})\in\tilde{\mathcal{F}}_{\mathcal{H}},$
jump to step $10,$ else set $k=k-1$ and go to step $8.$

Step 8: Calculate $a=\text{min}(\kappa_{1}^{k},\|\tilde{x}-\tilde{x}_{c}\|).$
If $a\leq\eta_{2},$ then go to step $9$ else jump back to step $5.$

Step 9: If $\tilde{u}_{1}\leq\epsilon_{2},$ then stop with a Probable Solution
else, set $i_{c}=i_{c}+1$ and jump back to Step $1$, after changing the
Initial Point as, $\tilde{x}^{(0)}=\tilde{x}_{c}.$

Step 10: Set $\tilde{t}_{s}=\tilde{\lambda}^{-1}(\tilde{u}_{s}),$
$\left[\begin{array}[]{c}\tilde{x}\\\ \tilde{t}\\\ \end{array}\right]=$
$\left[\begin{array}[]{c}\tilde{x}_{c}\\\ \tilde{t}_{s}\\\
\end{array}\right],\tilde{u}_{1}=\frac{\|\tilde{u}_{s}\|}{\sqrt{n}}$ If
$\tilde{u}_{1}\leq\epsilon_{1},$ then stop with acceptable homotopy solution
else set $i=i+1$ and go to step $2.$

######  Theorem 3.6:

If the homotopy curve $\Gamma_{\tilde{x}}^{(0)}$ is smooth, then the positive
predictor direction $\tilde{\tau}^{(0)}$ at the initial point
$\tilde{x}^{(0)}$ satisfies
sign($\det\left[\begin{array}[]{c}\frac{\partial\mathcal{H}}{\partial\tilde{x}\partial\tilde{\lambda}}(\tilde{x}^{(0)},1)\\\
e\tilde{\tau}^{(0)^{t}}\\\ \end{array}\right]$)$=(-1)^{n}\text{sign}\det
J_{\mathcal{H}}{\tilde{x}^{(0)}},$ where
$J_{\mathcal{H}}{\tilde{x}^{(0)}}=\frac{\partial{\psi}(\tilde{x}^{(0)})}{\partial\tilde{x}^{(0)}}.$

###### Proof.

From equation 3.4 we have
$\mathcal{H}(\tilde{x},\tilde{x}^{(0)},\tilde{\lambda})=$
$\left[\begin{array}[]{c}\psi_{1}(\tilde{x})-{\lambda_{1}}\psi_{1}(\tilde{x}^{(0)})\\\
\psi_{2}(\tilde{x})-{\lambda_{2}}\psi_{2}(\tilde{x}^{(0)})\\\ \vdots\\\
\psi_{n}(\tilde{x})-{\lambda_{n}}\psi_{n}(\tilde{x}^{(0)})\end{array}\right]=0.$
Now
$\frac{\partial\mathcal{H}}{\partial\tilde{x}\partial\tilde{\lambda}}(\tilde{x}=\tilde{x}^{(0)},\tilde{\lambda}=e)=$$\begin{bmatrix}J_{\mathcal{H}}\tilde{x}^{(0)}&\Psi(\tilde{x}^{(0)})\\\
\end{bmatrix},$
where
$\Psi(\tilde{x}^{(0)})=\begin{bmatrix}-\psi_{1}(\tilde{x}^{(0)})&0&0&\cdots&0\\\
0&-\psi_{2}(\tilde{x}^{(0)})&0&\cdots&0\\\
\vdots&\vdots&\vdots&\vdots&\vdots\\\
0&0&0&\cdots&-\psi_{n}(\tilde{x}^{(0)})\\\ \end{bmatrix}.$ Let positive
predictor direction be $\tau^{(0)}=\left[\begin{array}[]{c}\tilde{t}\\\
-e\end{array}\right]=\left[\begin{array}[]{c}(\tilde{R}^{(0)}_{1})^{(-1)}\tilde{R}_{2}^{(0)}e\\\
-e\end{array}\right],$ where
$\tilde{R}^{(0)}_{1}=J_{\mathcal{H}}\tilde{x}^{(0)}$ and
$\tilde{R}^{(0)}_{2}=\Psi(\tilde{x}^{(0)}).$ Here
$\text{det}(\tilde{R}^{(0)}_{1})\neq 0.$ Therefore,
$\det\left[\begin{array}[]{c}\frac{\partial\mathcal{H}}{\partial\tilde{x}^{(0)}\partial\tilde{\lambda}}\\\
e\tilde{\tau}^{(0)^{t}}\\\
\end{array}\right]$$=\det\left[\begin{array}[]{cc}R^{(0)}_{1}&R^{(0)}_{2}\\\
ee^{t}(R^{(0)}_{2})^{t}(R^{(0)}_{1})^{(-t)}&-ee^{t}\\\ \end{array}\right]$
$=\det\left[\begin{array}[]{cc}R^{(0)}_{1}&R^{(0)}_{2}\\\
0&-ee^{t}-ee^{t}(R^{(0)}_{2})^{t}(R^{(0)}_{1})^{(-t)}(R^{(0)}_{1})^{(-1)}R_{2}^{(0)}\\\
\end{array}\right]\\\
=\det(R^{(0)}_{1})\det(-ee^{t}-ee^{t}(R^{(0)}_{2})^{t}(R^{(0)}_{1})^{(-t)}(R^{(0)}_{1})^{(-1)}R_{2}^{(0)})\\\
=(-1)^{n}\det(R^{(0)}_{1})\det(ee^{t}(I+\mathcal{A}^{t}\mathcal{A})),$ where
$\mathcal{A}=(R^{(0)}_{1})^{(-1)}R_{2}^{(0)}.$ So
sign($\det\left[\begin{array}[]{c}\frac{\partial\mathcal{H}}{\partial\tilde{x}\partial\tilde{\lambda}}(\tilde{x}^{(0)},1)\\\
e\tilde{\tau}^{(0)^{t}}\\\ \end{array}\right]$)$=(-1)^{n}\text{sign}\det
J_{\mathcal{H}}{\tilde{x}^{(0)}}.$

##  4 Numerical Example

To illustrate the use of the above three algorithms presented in the previous
section, a numerical example[16] is presented here.
Following Murphy et al. [16] the oligopolistic market equilibrium problem
follows:

Let there be $n$ firms, which supply a homogeneous product in a noncooperative
fashion. Let $P(\tilde{Q}),\tilde{Q}\geq 0$ denote the inverse demand, where
$\tilde{Q}=\sum_{i=1}^{n}Q_{i},$ $Q_{i}\geq 0$ denote the $i$th firm’s supply.
Let $c_{i}(Q_{i})$ be the total cost of supplying $Q_{i}$ units. Now the Nash
equilibrium solution is a set of nonnegative output levels
${Q_{1}}^{*},{Q_{2}}^{*},\cdots,{Q_{n}}^{*},$ such that ${Q_{i}}^{*}$ is an
optimal solution to the following problem $\forall i\in\\{1,2\cdots,n\\}:$

$\displaystyle\text{maximize}\
{Q_{i}P(Q_{i}+{\tilde{Q}_{i}}^{*})-c_{i}(Q_{i})}$ (4.9)

where ${\tilde{Q}_{i}}^{*}=\sum_{j\neq i}{Q_{j}}^{*}.$ Murphy et al. show that
if $c_{i}(Q_{i})$ is convex and continuously differentiable $\forall
i\in\\{1,2,\cdots,n\\}$ and the inverse demand function $P(\tilde{Q})$ is
strictly decreasing and continuously differentiable and the industry revenue
curve $\tilde{Q}P(\tilde{Q})$ is concave, then
$({Q_{1}}^{*},{Q_{2}}^{*},\cdots,{Q_{n}}^{*})$ is a Nash equilibrium solution
if and only if

$\displaystyle[P(\tilde{Q}^{*})+{Q_{i}}^{*}P^{\prime}(\tilde{Q}^{*})-{c_{i}}^{\prime}({Q_{i}}^{*})]{Q_{i}}^{*}=0$
(4.10)
$\displaystyle{c_{i}}^{\prime}({Q_{i}}^{*})-P(\tilde{Q}^{*})-{Q_{i}}^{*}P^{\prime}(\tilde{Q}^{*})\geq
0$ (4.11) $\displaystyle{Q_{i}}^{*}\geq 0\ \ \forall i\in\\{1,2,\cdots,n\\}$
(4.12)

where $\tilde{Q}^{*}=\sum_{i=1}^{n}{Q_{i}}^{*},$ which is a nonlinear
complementarity problem with
$f_{i}(z)={c_{i}}^{\prime}({Q_{i}}^{*})-P(\tilde{Q}^{*})-{Q_{i}}^{*}P^{\prime}(\tilde{Q}^{*}),$
and $z_{i}={Q_{i}}^{*}.$
Note that here the functions $c_{i}(Q_{i})$ and $-\tilde{Q}P(\tilde{Q})$ are
convex. So the 1st order derivative of these two functions are increasing
function. Hence the function
$f_{i}(z)={c_{i}}^{\prime}({Q_{i}}^{*})-P(\tilde{Q}^{*})-{Q_{i}}^{*}P^{\prime}(\tilde{Q}^{*})$
is an increasing function.

Now consider an oligopoly with five firms, each with a total cost function of
the form:

$c_{i}(Q_{i})=n_{i}Q_{i}+\frac{\beta_{i}}{\beta_{i}+1}{L_{i}}^{\frac{1}{\beta_{i}}}{Q_{i}}^{\frac{\beta_{i}+1}{\beta_{i}}}$
(4.13)

The demand curve is given by:

$\tilde{Q}=5000P^{-1.1},\ \ P(\tilde{Q})=5000^{1/1.1}\tilde{Q}^{-1/1.1}.$
(4.14)

The parameters of the equation 4.13 for the five firms are given below:

Table 1: Value of parameters for five firms firm $i$ | $n_{i}$ | $L_{i}$ | $\beta_{i}$
---|---|---|---
1 | 10 | 5 | 1.2
2 | 8 | 5 | 1.1
3 | 6 | 5 | 1
4 | 4 | 5 | 0.8
5 | 2 | 5 | 0.6

To solve this problem using the homotopy method 3.4 with vector parameter
$\lambda,$ we first take the initial point
$\tilde{x}^{(0)}=$$\left[\begin{array}[]{c}1\\\ 1\\\ 1\\\ 1\\\ 1\\\
\end{array}\right],$ $\tilde{\lambda}^{(0)}=$$\left[\begin{array}[]{c}1\\\
1\\\ 1\\\ 1\\\ 1\\\ \end{array}\right].$ After $20$ iterations we obtain the
result $\tilde{x}=$ $\left[\begin{array}[]{c}15.429308\\\ 12.498582\\\
9.663473\\\ 7.165093\\\ 5.132566\\\ \end{array}\right],$
$\tilde{\lambda}=$$\left[\begin{array}[]{c}0\\\ 0\\\ 0\\\ 0\\\ 0\\\
\end{array}\right].$

##  5 Conclusion

In this study we show the equivalency between the nonlinear complementarity
problem and the system of nonlinear equations. We introduce a homotopy with
vector parameter to find the solution of the system of nonlinear equations as
well as the solution of the corresponding nonlinear complementarity problem.
In this connection we obtain the sign of the positive tangent direction of the
homotopy path. We construct a homotopy continuation method and show that the
method reaches to solution through a bounded and smooth curve under the
condition of the nonlinear function corresponding to the nonlinear
complementarity problem being either increasing function or bounded function.
Finally a real life oligopolistic market equilibrium problem is considered to
illustrate our results.

##  6 Acknowledgment

The author A. Dutta is thankful to the Department of Science and Technology,
Govt. of India, INSPIRE Fellowship Scheme for financial support.

## References

  * [1] Shui-Nee Chow, John Mallet-Paret, and J. A. Yorke.  Finding zeroes of maps : Homotopy methods that are constructive with probability one *.  2010\.
  * [2] AK Das.  Properties of some matrix classes based on principal pivot transform.  Annals of Operations Research, 243, 05 2014.
  * [3] AK Das, R Jana, and Deepmala.  On generalized positive subdefinite matrices and interior point algorithm.  In Frontiers in Optimization: Theory and Applications, pages 3–16. Springer, 2016.
  * [4] AK Das, R Jana, and Deepmala.  Finiteness of criss-cross method in complementarity problem.  In International Conference on Mathematics and Computing, pages 170–180. Springer, 2017.
  * [5] A Dutta, R Jana, and AK Das.  On column competent matrices and linear complementarity problem.  In Proceedings of the Seventh International Conference on Mathematics and Computing, pages 615–625. Springer, 2022.
  * [6] B Curtis Eaves and Romesh Saigal.  Homotopies for computation of fixed points on unbounded regions.  Mathematical Programming, 3(1):225–237, 1972.
  * [7] Jin-yan Fan.  A modified levenberg-marquardt algorithm for singular system of nonlinear equations.  Journal of Computational Mathematics, pages 625–636, 2003.
  * [8] Hossein Jafari and Varsha Daftardar-Gejji.  Revised adomian decomposition method for solving a system of nonlinear equations.  Applied Mathematics and Computation, 175(1):1–7, 2006.
  * [9] R Jana, AK Das, and A. Dutta.  On hidden z-matrix and interior point algorithm.  OPSEARCH, 56, 09 2019.
  * [10] R Jana, A. Dutta, and AK Das.  More on hidden z -matrices and linear complementarity problem.  Linear and Multilinear Algebra, 69:1–10, 06 2019.
  * [11] Stepan Karamardian.  The nonlinear complementarity problem with applications, part 1.  Journal of Optimization Theory and Applications, 4(2):87–98, 1969\.
  * [12] SR Mohan, S.K. Neogy, and AK Das.  More on positive subdefinite matrices and the linear complementarity problem.  Linear Algebra and Its Applications, 338(1-3):275–285, 2001.
  * [13] SR Mohan, S.K. Neogy, and AK Das.  On the classes of fully copositive and fully semimonotone matrices.  Linear Algebra and its Applications, 323:87–97, 01 2001.
  * [14] S.R. Mohan, S.K. Neogy, and AK Das.  A note on linear complementarity problems and multiple objective programming.  Mathematical Programming. Series A. Series B, 100, 06 2004.
  * [15] Prasenjit Mondal, S Sinha, Samir K Neogy, and AK Das.  On discounted ar–at semi-markov games and its complementarity formulations.  International Journal of Game Theory, 45(3):567–583, 2016.
  * [16] Frederic H Murphy, Hanif D Sherali, and Allen L Soyster.  A mathematical programming approach for determining oligopolistic market equilibrium.  Mathematical Programming, 24(1):92–106, 1982.
  * [17] Samir K Neogy, AK Das, S Sinha, and A Gupta.  On a mixture class of stochastic game with ordered field property.  In Mathematical programming and game theory for decision making, pages 451–477. World Scientific, 2008.
  * [18] S.K. Neogy, R Bapat, AK Das, and T Parthasarathy.  Mathematical Programming and Game Theory for Decision Making.  04 2008.
  * [19] S.K. Neogy, R. Bapat, AK Das, and T. Parthasarathy.  Mathematical programming and game theory for decision making.  11 2021.
  * [20] S.K. Neogy, Ravindra Bapat, AK Das, and Biswabrata Pradhan.  Optimization models with economic and game theoretic applications.  Annals of Operations Research, 243, 07 2016.
  * [21] S.K. Neogy and AK Das.  Linear complementarity and two classes of structured stochastic games.  Operations Research with Economic and Industrial Applications: Emerging Trends, eds: SR Mohan and SK Neogy, Anamaya Publishers, New Delhi, India, pages 156–180, 2005.
  * [22] S.K. Neogy and AK Das.  On almost type classes of matrices with q-property.  Linear & Multilinear Algebra - LINEAR MULTILINEAR ALGEBRA, 53:243–257, 07 2005.
  * [23] S.K. Neogy and AK Das.  Some properties of generalized positive subdefinite matrices.  SIAM J. Matrix Analysis Applications, 27:988–995, 01 2006.
  * [24] S.K. Neogy and AK Das.  On weak generalized positive subdefinite matrices and the linear complementarity problem.  Linear and Multilinear Algebra, 61, 07 2013.
  * [25] S.K. Neogy, AK Das, and R. Bapat.  Modeling, computation and optimization.  06 2022.
  * [26] S.K. Neogy, AK Das, and Abhijit Gupta.  Generalized principal pivot transforms, complementarity theory and their applications in stochastic games.  Optimization Letters, 6:339–356, 02 2012.
  * [27] S.K. Neogy, AK Das, and Abhijit Gupta.  Generalized principal pivot transforms, complementarity theory and their applications in stochastic games.  Optimization Letters, 6(2):339–356, 2012.
  * [28] Muhammad Aslam Noor and Muhammad Waseem.  Some iterative methods for solving a system of nonlinear equations.  Computers & Mathematics with Applications, 57(1):101–106, 2009\.
  * [29] Takeo Ojika, Satoshi Watanabe, and Taketomo Mitsui.  Deflation algorithm for the multiple roots of a system of nonlinear equations.  Journal of mathematical analysis and applications, 96(2):463–479, 1983.
  * [30] K Sayevand and H Jafari.  On systems of nonlinear equations: some modified iteration formulas by the homotopy perturbation method with accelerated fourth-and fifth-order convergence.  Applied Mathematical Modelling, 40(2):1467–1476, 2016.
  * [31] Xiuyu Wang and Xingwu Jiang.  A homotopy method for solving the horizontal linear complementarity problem.  Computational and Applied Mathematics, 33, 04 2013.
  * [32] Layne T Watson.  Globally convergent homotopy methods: a tutorial.  Applied Mathematics and Computation, 31:369–396, 1989.
  * [33] Layne T Watson and Raphael T Haftka.  Modern homotopy methods in optimization.  Computer Methods in Applied Mechanics and Engineering, 74(3):289–305, 1989.
  * [34] Kiyotaka Yamamura, Hitomi Kawata, and Ai Tokue.  Interval solution of nonlinear equations using linear programming.  BIT Numerical Mathematics, 38(1):186–199, 1998.
  * [35] X Zhao, S Zhang, and Q Liu.  A combined homotopy interior point method for the linear complementarity problem.  Journal of Information and Computational Science, 7:1589–1594, 07 2010.