11institutetext: School of Engineering, The University of Edinburgh,
Edinburgh, EH9 3FG, UK 22institutetext: Krannert Cardiovascular Research
Center, Indiana University School of Medicine, Indianapolis, Indiana,
USA33institutetext: The Alan Turing Institute, London, NW1 2DB, UK

# The MRI Scanner as a Diagnostic:
Image-less Active Sampling

Yuning Du 11 Rohan Dharmakumar 22 Sotirios A.Tsaftaris 11 3 3

###### Abstract

Despite the high diagnostic accuracy of Magnetic Resonance Imaging (MRI),
using MRI as a Point-of-Care (POC) disease identification tool poses
significant accessibility challenges due to the use of high magnetic field
strength and lengthy acquisition times. We ask a simple question: Can we
dynamically optimise acquired samples, at the patient level, according to an
(automated) downstream decision task, while discounting image reconstruction?
We propose an ML-based framework that learns an active sampling strategy, via
reinforcement learning, at a patient-level to directly infer disease from
undersampled $k$-space. We validate our approach by inferring Meniscus Tear in
undersampled knee MRI data, where we achieve diagnostic performance comparable
with ML-based diagnosis, using fully sampled $k$-space data. We analyse task-
specific sampling policies, showcasing the adaptability of our active sampling
approach. The introduced frugal sampling strategies have the potential to
reduce high field strength requirements that in turn strengthen the viability
of MRI-based POC disease identification and associated preliminary screening
tools.

###### Keywords:

Point-of-Care Diagnosis Active Sampling StrategyReinforcement Learning
Magnetic Resonance Imaging.

##  1 Introduction

Despite the proliferation of Magnetic Resonance Imaging (MRI), its role as a
Point-Of-Care (POC) diagnostic tool is muted due to poor accessibility caused
by long acquisition time and cumbersome equipment [4]. Consequently,
advancements have been made in assisting the diagnostic process with machine
learning methods to sample less $k$-space data and reduce acquisition time in
turn [8].

We first illustrate how learning-based strategies can help to reduce
acquisition time by considering the depiction of conventional MRI-based
diagnostic processes. In Figure 1a. an MRI scanner samples the full $k$-space,
resulting in high-fidelity images, which are further interpreted by
professional radiologists to identify biomarkers and provide a diagnosis.
Alternatively, machine-learning (ML) based diagnostic processes can enable
acquisition time savings by reconstructing a high-fidelity image from
undersampled $k$-space (Figure 1b). Previous work has focused on optimising
both reconstruction [3] and $k$-space sampling strategies [1]. To find the
best $k$-space sampling pattern, various methods [17, 18, 12] optimise the
mask given pre-defined sample rates, resulting in population-level masks,
while reinforcement-learning-based methods [2, 11] can optimise an active
sampling strategy at population or patient-level. A re-interpretation of such
conventional use of MRI can lead to considerable savings, which will largely
enable a low-field MRI future opening the road to POC and bedside imaging.

Figure 1: Conventional vs. ML-based Diagnostic Processes.

Considering that image reconstruction is a shared stage in both conventional
(Figure 1a) and ML-based (Figure 1b) diagnostic processes, one ponders Is the
image necessary for automated, diagnostic, inferences? One option is to
perform direct inference on undersampled $k$-space as shown in Figure 1c. In
fact, [14] have shown that it is possible to obtain biomarkers directly from
$k$-space data. However, one then immediately ponders Can AI agents learn
effective $k$-space sampling strategies by considering expected diagnostic
performance?  A recent realisation of this concept is presented in [15]. The
work formulates the described process as a classification task with an
optimised sampling strategy. However, the sampling strategy is designed at the
population-level, resulting in the same mask for all patients. All approaches
that find optimal patient-level active sampling strategies use image
reconstruction [2, 11, 20]. To the best of our knowledge, active sampling
strategies for patient-level disease inference, from undersampled $k$-space,
are lacking. Motivated by this, our contributions include:

  1. $\bullet$

We propose, the first, patient-level active sampling using a reinforcement
learning framework, aiming for direct disease inference from $k$-space.

  2. $\bullet$

We test feasibility to infer Meniscus Tear presence in undersampled MRI.

  3. $\bullet$

We investigate how different policies make decisions and how different
starting points alter behaviour.

##  2 MRI Undersampling Preliminaries

Instead of directly imaging the human anatomy, MRI captures the
electromagnetic activity in the body after exposure to magnetic fields and
radiofrequency pulses, which can be measured in $k$-space (i.e., the frequency
domain). Considering the single coil measurement, the $k$-space data can be
represented as a 2-dimensional complex-valued matrix $x\in\mathbb{C}^{r\times
c}$, where $r$ is the number of rows and $c$ is the number of columns. The
spatial image $I$ can be obtained by applying the inverse Fourier Transform to
$x$, denoted as $I=\mathcal{F}^{-1}(x)$. The undersampled $k$-space can be
therefore represented as $x_{s}=U_{L}\circ x$, where $U_{L}$ can be viewed as
a binary mask $U\in\\{0,1\\}^{r\times c}$ with L measurements from $k$-space.
In our work, we exclusively consider the Cartesian mask for MRI undersampling.
Consequently, the undersampled image is denoted as
$I_{s}=\mathcal{F}^{-1}(x_{s})$.

Figure 2: Our framework for active MRI sampling for point-of-care diagnosis.

##  3 Methods

Overview: Our framework, illustrated in Figure 2, aims to reduce acquisition
time by selectively sampling the $k$-space according to diagnostic
significance, in a progressive fashion. The process begins with a small subset
of randomly sampled $k$-space with L measurements denoted as
$x_{s_{0}}=U_{L}\circ x$, resulting in a low-quality image $I_{s_{0}}$.
$I_{s_{0}}$ is then input into a pre-trained classification network $f$ to
generate the initial prediction $y_{0}$. The extracted feature maps $m_{0}$
are processed by the Feature Operator $O$, producing high-level features
$h_{0}=O(m_{0})$, which are then fed into the Active Sampler $S$. The policy
network generates a sampling policy $p_{\phi}(h_{0})$ parameterised by $\phi$,
guiding the selective sampling of diagnostically significant lines from the
$k$-space. These sampled lines are subsequently added to $U_{L}$. The updated
undersampled $k$-space data at step $t\in[0,T]$ are denoted as
$x_{s_{t}}=U_{t}\circ x$, and are then fed back into the classification
network with inverse Fourier transform. The iterative process continues until
a sampling budget $T$ is exhausted or user-defined reliability criteria are
satisfied. During the training stage, the predictions $y_{t}$, accompanied by
the ground truth diagnostic label $g$, are used as the criterion
$r_{t}=r(y_{t},g)$ to supervise the $S$. While in inferece stage, there is no
ground truth diagnostic label $g$ provided.

Classification Network and Image Operator: To improve stability, the
classification network is pre-trained with undersampled $k$-space data,
ensuring an advantageous reward for the active sampler during training. This
network also functions as the feature extractor in our setting. In the
classification network, the earlier layers tend to learn low-level features
such as edges, textures, and simple patterns, while deeper layers learn more
complex and high-level features that are useful for discriminating between
different classes. Therefore, utilising the feature maps from the model, the
high-level features are selected and processed by the feature operator as
further inputs.

Active $k$-space Sampler with Greedy Policy: Inspired by [2, 11], the
sequential selection of $k$-space can be formalised as a Partially Observable
Markov Decision Process (POMDP) [16]. Greedy Policy is used to maximise the
expected return $J(\phi)$ of a policy $p_{\phi}$ parameterised by $\phi$ in
such a POMDP. At each step, the classification improvement can be calculated
using $R(r_{t},r_{t+1})=r(y_{t+1},g)-r(y_{t},g)$, where the criterion $r$ is
the cross-entropy.

During inference, the agent will be sampling one-line at a time. However, this
will slow during training. Hence, the policy network is trained by sampling
several lines in parallel the rewards of which are averaged [6]. Formally, we
sample $q$ lines at every time step, for a reward $R_{i,t}$ as the reward
obtained from sample $i$ at time step $t$, to obtain the following estimator:

$\nabla_{\phi}J(\phi)\approx\frac{1}{q-1}\mathbb{E}_{x}\sum_{i=1}^{q}\sum_{t=L}^{T-1}\left[\nabla_{\phi}\log
p_{\phi}\left(h_{t}\right)\left(R_{i,t}-\frac{1}{q}\sum_{j=1}^{q}R_{j,t}\right)\right].$

(1)

Evaluation Metrics: To assess the network’s classification performance, we
employ metrics such as Recall, Area Under Curve (AUC), and Specificity.

##  4 Experiments

###  4.1 Dataset and Pre-processing

Dataset: We used single-coil $k$-space data and slice-level labels from the
publicly available fastMRI dataset [19] and fastMRI+ dataset [21]. Randomly
selecting $1100$ annotated volumes ($39125$ slices) from the fastMRI Knee
dataset. Our diagnostic task is to identification Meniscus Tear (MT) in each
slice. Thereby, there are $32035$ train slices ($11.8\%$ with MT), $5249$
validation slices ($9.3\%$ with MT), and $1841$ test slices ($11.1\%$ with
MT).

Table 1: Diagnosis Support and Data Access of Implemented Methods.

| Diagnostic Support | Data Used by the Model
---|---|---
Method |  Sampling Optimization |  Patient-level Strategy |  Full $k$-space |  Diagnostic Label |  Undersampled $k$-space
Oracle |  |  | ✓ | ✓ |
Undersampled |  |  |  | ✓ | ✓
Policy Reconstruction [2] | ✓ | ✓ | ✓ | ✓ | ✓
Policy Classifier(Ours) | ✓ | ✓ |  | ✓ | ✓

Data Pre-processing: Since the $k$-space data have various sizes, we first use
inverse Fourier transform to the fully sampled $k$-space data to get the
ground truth image and crop it to size $(320\times 320)$ for computation
convenience. Thereby, the fully sampled $k$-space data of uniform size can be
obtained by applying Fourier transform to the ground truth image. Notably,
there is a severe class imbalance regarding the MT identification task. During
training the classification network, we oversample the data to avoid
overfitting on the majority class and poor generalisation on the minority
class.

###  4.2 Implementation Details

We devise inference benchmarks to allow us to evaluate our approach in a fair
fashion. Two benchmarks have the access to fully sampled $k$-space(high-
fidelity images), and two do not optimise at the patient level the sampling
strategy. The diagnosis support and data access are shown in Table 1.

Fully Sampled (Oracle): This serves as an benchmark estimator of classifier
performance on image input obtained by fully sampled $k$-space data, and hence
no sampling optimisation occurs. We trained the classifier with the ground
truth image as input which is transformed from fully sampled $k$-space data
and supervised by the diagnostic label using cross-entropy loss. We use as
classification backbone a ResNet-50 [5] and to address the class imbalance in
the training set, we add extra dropout layers overfitting, resulting in a
total of $25.6\mathrm{M}$ parameters.

Undersampled: This classifier serves as a baseline of performance when simple
inverse Fourier is used to transform the under-sampled $k$-space data, without
any sampling optimisation. It has the same backbone as the Oracle, and is
trained with undersampled images with various sample rates ($4$ to $20$) and
center fraction ($0$ to $0.10$), supervised as before.

Policy (via) Reconstruction [2]: We compare with a model that optimises
patient-level sampling strategy with image reconstruction error as rewards.
Notably this method has access to fully sampled $k$-space, and hence has
access to more information during training. We pre-trained a reconstruction
network using a U-Net [13] as the backbone with a first feature map size of
$16$ and $4$ pooling cascades, resulting in a total of $837\mathrm{K}$
parameters. The reconstruction network is trained with various sample rates
($4$ to $20$) and center fraction ($0$ to $0.10$) supervised using the
$\ell_{1}$ loss. We train the active sampler for reconstruction with
reconstructed images from pre-trained model as inputs, and Structural
Similarity Index Measure as the criterion to provide rewards, resulting in a
total of $26.7\mathrm{M}$ parameters. The reconstructed images obtained with
the policy network and pre-trained reconstruction model are evaluated with the
Oracle, referred to as the Policy-based Reconstruction Network (Policy
Reconstruction).

Proposed Policy Classifier: For our method, a backbone similar to the
Undersampled is used as pre-trained classification network. The reward is
driven by the predictions of the network and its feature maps are used to
train the policy. The Feature Operator uses the last 2 layers’ feature maps as
input and the global average pooling function to achieve a $(80,512)$ output
to feed in the policy network. The policy network ($11\mathrm{M}$ parameters)
uses the cross entropy from classification network to provide rewards. It is
trained with an initial sample rate of $5\%$ with multiple center fraction and
samples $64$ lines to reach the sample rate of $25\%$ as the sampling budget.
The parallel acquisition $q=8$.

For all methods, we employ the Adam optimiser with a learning rate of
$10^{-4}$ and a step-based scheduler with a decay gamma of $0.1$ for all model
training. All classification and reconstruction models are trained for $30$
epochs, and the policy networks of the active sampler are trained for $20$
epochs. Our experimental setup uses the PyTorch framework, and all
computations are conducted on NVIDIA A100 Tensor Core GPUs. Our
code111https://anonymous.4open.science/r/KspaceToDiagnosis-CB5E is available.

###  4.3 Results

Figure 3: Classification performance varying sample rates (horizontal axis).
More metrics are reported in the Supplementary.

Making Diagnostic Decisions with Undersampled Data
Figure 3 compares the performance at various sample rates for three strategies
that use undersampled masks, namely the Undersampled, Policy Reconstruction
and Policy Classifier. The first randomly samples $k$-space lines
progressively; the other two start with a randomly initialised mask with a
sample rate of $5\%$ and optimally decide using their respective rewards and
policies (see Section 3). Our method, the Policy Classifier, consistently
outperforms the Undersampled baseline across all metrics, showing that optimal
sampling via the policy network helps in identifying diagnostically relevant
$k$-space lines. This leads to considerable savings in data sampled (and
consequently scan time). Our model approximates well the performance of the
Oracle, which has been trained on high-fidelity data, and reaches optimal
performance quickly. Taking AUC as an example, our approach reaches an AUC of
0.780 with 7% of the sampled $k$-space. Our method closely approximates the
performance of the Policy Reconstruction, which we highlight has been trained
with access to fully sampled $k$-space.222 We observe that the Policy
Reconstruction performs better than the Oracle. This policy reconstructs an
image which is given to the Oracle for classification. Some denoising and
smoothing are happening at the reconstruction which in turn acts as a
regulariser for the Oracle classifier explaining this slightly improved
performance.

Figure 4: $k$-space behaviour of the two policies. The horizontal axis
indicates the cumulative lines acquired while sampling.The initial 16 lines
are randomly sampled.

$k$-space Sampling Behaviour of the Policies
The results of the previous paragraph are reported over a coarse percentage of
$k$-space lines acquired. It is worth looking into how the two different
policies behave when asking each policy to progressively make decisions on
which line to acquire in a line-by-line manner. The results of this exercise
are shown in Figure 4. The behaviour of the two methods is evidently
different. Taking Recall and AUC as examples, the Policy Classifier plateaus
quickly, reaching excellent performance with 22 sampled lines and continuing
to make, small, improvements. The Policy Reconstruction over the same interval
makes sub-optimal (to classification) decisions before it plateaus. This
behaviour can be explained by the fact that a reconstruction policy may not
favour lines useful for classification. The Policy Classifier appears to
change behaviour between 19 and 21 lines acquired. This can be attributed to a
drift from how the classifier was pre-trained with randomly sampled lines. We
discuss solutions to this later.

Figure 5: $k$-space preference as a function of sample rate (vertical axis
downward) and for 3 center fraction scenarios. Colorbar indicates the
possibility of being sampled.

Task-specific ‘Coarse-to-Fine’ Sampling Policy
It would be interesting to see which $k$-space lines are favoured by our
policy and how altering the starting point (namely different percentage of
center fraction which represents the amount of low-frequency $k$-space
intentionally sampled from the center of the $k$-space) modulates this
behaviour. Figure 5 illustrates the average masks provided by three policy
networks. Comprared with the sampling policy for reconstruction in
supplementary, this analysis reveals that early on the policy favours both
low-frequency and some high-frequency data. However, later on the policy
favours specific diagnostically-relevant $k$-space lines of high-frequency.
Hence we see the policy capturing ‘Coarse’ features such as anatomical
structures or essential patterns first and later ‘Fine’ details from the high-
frequency data. When center fraction increases, the policy samples even more
high-frequency information early on. While this might be obvious it is
actually driven by low AUC performance when forcing the model to sample more
center lines from the start (the AUC drops from $0.794$ to $0.743$ when the
center fraction is set to $0.05$, see also more metrics in the supplementary).
By forcing the model to sample redundant low-frequency data, the policy tries
to recover to identify diagnostically significant data within the $k$-space.

##  5 Discussion

Clinical Relevance: Our method actively samples $k$-space data whilst
simultaneously conducting disease inference, enabling real-time diagnosis
during scanning. This approach eliminates the need for sampling full
$k$-space, thereby reducing acquisition time. In addition, our approach does
not require high-fidelity fully sampled data for policy training, which is a
hard requirement of reconstruction policies. We believe this has the potential
to further accelerate the development of low-field MR. We envision
applications in pre-screening by offering preliminary diagnostic results to
aid in resource allocation, particularly in regions with limited medical
resources.

Limitations: Our proof of concept is based on the fastMRI knee data and
specifically MT. One could readily see the application of a portable low-field
scanner in sports medicine [10] but it is a less compelling clinical
application. Ideally, a dataset of e.g. stroke and acute brain trauma [9],
which would probe diverse, and not only structural, MR signal characteristics
would be more compelling. At this moment such dataset is currently lacking.
Our policy network relies on a pre-trained classifier, which requires pre-
training with diverse undersampled data. It is of great interest to train the
classifier simultaneously with the policy network. Finally, our policy network
does not leverage physical priors of $k$-space utility which are implicitly
leveraged by reconstruction policies using a pre-trained reconstruction
network. Such prior also may encode physical limitations [7] of the employed
sequence in making sampling decisions.

##  6 Conclusion

Our proposed framework presents a novel approach to enhance the efficiency and
accessibility of MRI as a Point of Care (POC) diagnostic tool by optimising
sampling patterns in a learnable fashion. Our approach distinctly stipulates
that $k$-space sampling acquisition decisions are undertaken by an agent
optimising a dynamic per-patient classification task. Results reporting direct
inference from undersampled $k$-space data, concerning the presence of
Meniscus Tear, showed that our approach achieves the comparable performance as
an Oracle Classifier whilst reducing acquisition time with only 25% $k$-space
usage. Furthermore, our analysis of task-specific sampling policies revealed
that the policy network adapted its sampling strategy based on the nature of
the task, demonstrating the adaptability of the active sampling approach.

## References

  * [1] Bahadir, C.D., Wang, A.Q., Dalca, A.V., Sabuncu, M.R.: Deep-learning-based optimization of the under-sampling pattern in mri. IEEE Transactions on Computational Imaging 6, 1139–1152 (2020)
  * [2] Bakker, T., van Hoof, H., Welling, M.: Experimental design for mri by greedy policy search. Advances in Neural Information Processing Systems 33, 18954–18966 (2020)
  * [3] Cai, L., Gao, J., Zhao, D.: A review of the application of deep learning in medical image classification and segmentation. Annals of translational medicine 8(11) (2020)
  * [4] Geethanath, S., Vaughan Jr, J.T.: Accessible magnetic resonance imaging: a review. Journal of Magnetic Resonance Imaging 49(7), e65–e77 (2019)
  * [5] He, K., Zhang, X., Ren, S., Sun, J.: Identity mappings in deep residual networks. In: Computer Vision–ECCV 2016: 14th European Conference, Amsterdam, The Netherlands, October 11–14, 2016, Proceedings, Part IV 14. pp. 630–645. Springer (2016)
  * [6] Kool, W., van Hoof, H., Welling, M.: Buy 4 reinforce samples, get a baseline for free! (2019)
  * [7] Kuperman, V.: Magnetic resonance imaging: physical principles and applications. Elsevier (2000)
  * [8] Lin, D.J., Johnson, P.M., Knoll, F., Lui, Y.W.: Artificial intelligence for mr image reconstruction: an overview for clinicians. Journal of Magnetic Resonance Imaging 53(4), 1015–1028 (2021)
  * [9] Lyu, M., Mei, L., Huang, S., Liu, S., Li, Y., Yang, K., Liu, Y., Dong, Y., Dong, L., Wu, E.X.: M4raw: A multi-contrast, multi-repetition, multi-channel mri k-space dataset for low-field mri research. Scientific Data 10(1), 264 (2023)
  * [10] Massimiliano, L., Giuseppe, G., Michela, B., Michele, A., Silvia, R., Alessio, P., Francesco, P., Lia, R., Alessandro, S., Federico, A.G., et al.: Role of low field mri in detecting knee lesions. Acta Bio Medica: Atenei Parmensis 90(Suppl 1), 116 (2019)
  * [11] Pineda, L., Basu, S., Romero, A., Calandra, R., Drozdzal, M.: Active mr k-space sampling with reinforcement learning. In: Medical Image Computing and Computer Assisted Intervention–MICCAI 2020: 23rd International Conference, Lima, Peru, October 4–8, 2020, Proceedings, Part II 23. pp. 23–33. Springer (2020)
  * [12] Ravula, S., Levac, B., Jalal, A., Tamir, J.I., Dimakis, A.G.: Optimizing sampling patterns for compressed sensing mri with diffusion generative models. arXiv preprint arXiv:2306.03284 (2023)
  * [13] Ronneberger, O., Fischer, P., Brox, T.: U-net: Convolutional networks for biomedical image segmentation. In: Medical Image Computing and Computer-Assisted Intervention–MICCAI 2015: 18th International Conference, Munich, Germany, October 5-9, 2015, Proceedings, Part III 18. pp. 234–241. Springer (2015)
  * [14] Schlemper, J., Oktay, O., Bai, W., Castro, D.C., Duan, J., Qin, C., Hajnal, J.V., Rueckert, D.: Cardiac mr segmentation from undersampled k-space using deep latent representation learning. In: Medical Image Computing and Computer Assisted Intervention–MICCAI 2018: 21st International Conference, Granada, Spain, September 16-20, 2018, Proceedings, Part I. pp. 259–267. Springer (2018)
  * [15] Singhal, R., Sudarshan, M., Mahishi, A., Kaushik, S., Ginocchio, L., Tong, A., Chandarana, H., Sodickson, D.K., Ranganath, R., Chopra, S.: On the feasibility of machine learning augmented magnetic resonance for point-of-care identification of disease. arXiv preprint arXiv:2301.11962 (2023)
  * [16] Sondik, E.J.: The optimal control of partially observable Markov processes. Stanford University (1971)
  * [17] Wang, J., Yang, Q., Yang, Q., Xu, L., Cai, C., Cai, S.: Joint optimization of cartesian sampling patterns and reconstruction for single-contrast and multi-contrast fast magnetic resonance imaging. Computer Methods and Programs in Biomedicine 226, 107150 (2022)
  * [18] Xuan, K., Sun, S., Xue, Z., Wang, Q., Liao, S.: Learning mri k-space subsampling pattern using progressive weight pruning. In: Medical Image Computing and Computer Assisted Intervention–MICCAI 2020: 23rd International Conference, Lima, Peru, October 4–8, 2020, Proceedings, Part II 23. pp. 178–187. Springer (2020)
  * [19] Zbontar, J., Knoll, F., Sriram, A., Murrell, T., Huang, Z., Muckley, M.J., Defazio, A., Stern, R., Johnson, P., Bruno, M., et al.: fastmri: An open dataset and benchmarks for accelerated mri. arXiv preprint arXiv:1811.08839 (2018)
  * [20] Zhang, Z., Romero, A., Muckley, M.J., Vincent, P., Yang, L., Drozdzal, M.: Reducing uncertainty in undersampled mri reconstruction with active acquisition. In: Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition. pp. 2049–2058 (2019)
  * [21] Zhao, R., Yaman, B., Zhang, Y., Stewart, R., Dixon, A., Knoll, F., Huang, Z., Lui, Y.W., Hansen, M.S., Lungren, M.P.: fastmri+: Clinical pathology annotations for knee and brain fully sampled multi-coil mri data. arXiv preprint arXiv:2109.03812 (2021)

## Supplementary

Table 2: Image Classification Performance with Various Sample Rates and Center
Fraction. Numbers in bold represent the upper bound.

Center Fraction | 0.00 | 0.01 | 0.05
---|---|---|---
AUC$\uparrow$
Sample Rate | 5% | 7% | 10% | 13% | 17% | 25% | 5% | 7% | 10% | 13% | 17% | 25% | 5% | 7% | 10% | 13% | 17% | 25%
Oracle | 0.819 | 0.819 | 0.819 | 0.819 | 0.819 | 0.819 |  | 0.819 | 0.819 | 0.819 | 0.819 | 0.819 | 0.819 | 0.819 | 0.819 | 0.819 | 0.819 | 0.819
Undersampled | 0.666 | 0.689 | 0.710 | 0.736 | 0.745 | 0.756 | 0.729 | 0.761 | 0.789 | 0.764 | 0.776 | 0.783 | 0.752 | 0.751 | 0.747 | 0.745 | 0.743 | 0.749
Recon Classifier | 0.586 | 0.612 | 0.648 | 0.647 | 0.652 | 0.687 | 0.782 | 0.786 | 0.774 | 0.773 | 0.752 | 0.778 | 0.798 | 0.813 | 0.829 | 0.835 | 0.838 | 0.836
Policy Reconstruction | 0.583 | 0.806 | 0.829 | 0.844 | 0.841 | 0.837 | 0.757 | 0.808 | 0.817 | 0.826 | 0.827 | 0.813 | 0.789 | 0.818 | 0.833 | 0.829 | 0.827 | 0.831
Policy Classifier | 0.666 | 0.777 | 0.775 | 0.781 | 0.785 | 0.791 | 0.717 | 0.789 | 0.785 | 0.784 | 0.775 | 0.794 | 0.752 | 0.755 | 0.739 | 0.742 | 0.741 | 0.743
RECALL$\uparrow$
Sample Rate | 5% | 7% | 10% | 13% | 17% | 25% | 5% | 7% | 10% | 13% | 17% | 25% | 5% | 7% | 10% | 13% | 17% | 25%
Oracle | 0.816 | 0.816 | 0.816 | 0.816 | 0.816 | 0.816 | 0.816 | 0.816 | 0.816 | 0.816 | 0.816 | 0.816 | 0.816 | 0.816 | 0.816 | 0.816 | 0.816 | 0.816
Undersampled | 0.592 | 0.636 | 0.655 | 0.699 | 0.699 | 0.694 | 0.592 | 0.655 | 0.728 | 0.680 | 0.704 | 0.723 | 0.636 | 0.636 | 0.626 | 0.621 | 0.621 | 0.641
Recon Classifier | 0.243 | 0.311 | 0.388 | 0.379 | 0.408 | 0.510 | 0.762 | 0.772 | 0.767 | 0.762 | 0.733 | 0.791 | 0.748 | 0.791 | 0.840 | 0.854 | 0.869 | 0.869
Policy Reconstruction | 0.228 | 0.850 | 0.893 | 0.917 | 0.908 | 0.888 | 0.699 | 0.835 | 0.859 | 0.884 | 0.888 | 0.854 | 0.728 | 0.816 | 0.859 | 0.854 | 0.845 | 0.854
Policy Classifier | 0.592 | 0.699 | 0.694 | 0.709 | 0.718 | 0.738 | 0.544 | 0.714 | 0.709 | 0.709 | 0.694 | 0.709 | 0.636 | 0.646 | 0.617 | 0.621 | 0.621 | 0.626
SPECIFICITY$\uparrow$
Sample Rate | 5% | 7% | 10% | 13% | 17% | 25% | 5% | 7% | 10% | 13% | 17% | 25% | 5% | 7% | 10% | 13% | 17% | 25%
Oracle | 0.822 | 0.822 | 0.822 | 0.822 | 0.822 | 0.822 | 0.822 | 0.822 | 0.822 | 0.822 | 0.822 | 0.822 | 0.822 | 0.822 | 0.822 | 0.822 | 0.822 | 0.822
Undersampled | 0.740 | 0.741 | 0.765 | 0.774 | 0.791 | 0.817 | 0.866 | 0.866 | 0.851 | 0.848 | 0.848 | 0.843 | 0.869 | 0.865 | 0.869 | 0.869 | 0.864 | 0.857
Recon Classifier | 0.928 | 0.914 | 0.907 | 0.916 | 0.895 | 0.864 | 0.801 | 0.799 | 0.781 | 0.783 | 0.771 | 0.764 | 0.849 | 0.834 | 0.818 | 0.815 | 0.807 | 0.804
Policy Reconstruction | 0.938 | 0.763 | 0.765 | 0.770 | 0.774 | 0.785 | 0.815 | 0.780 | 0.774 | 0.768 | 0.766 | 0.771 | 0.850 | 0.821 | 0.807 | 0.804 | 0.810 | 0.808
Policy Classifier | 0.740 | 0.854 | 0.856 | 0.853 | 0.851 | 0.843 | 0.890 | 0.864 | 0.861 | 0.860 | 0.856 | 0.852 | 0.869 | 0.864 | 0.862 | 0.862 | 0.861 | 0.861

Figure 6: $k$-space preference for reconstruction task as a function of sample
rate (vertical axis downward) and for 3 center fraction scenarios. Colorbar
indicates the possibility of being sampled.


# Tracing homotopy path for the solution of nonlinear complementarity Problem

A. Duttaa,1, 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

In this article, we consider nonlinear complementarity problem. We introduce a
new homotopy function for finding the solution of nonlinear complementarity
problem through the trajectory . We show that the homotopy path approaching
the solution is smooth and bounded . Numerical example of an oligopoly
equilibrium problem is illustrated to show the effectiveness of the proposed
algorithm.
Keywords: Nonlinear complementarity problem, homotopy continuation method,
Oligopoly equilibrium problem.

##  1 Introduction

The nonlinear complementarity problem is identified as an important
mathematical programming problem. The idea of nonlinear complementarity
problem is based on the concept of linear complementarity problem. For recent
study on this problem and applications see [5], [27], [20], [22] and
references therein. For details of several matrix classes in complementarity
theory, see [10], [11], [26], [24], [17], [16], [21], [25], [6] 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 [18], [19],
[23] and [30]. The complementarity problems are considered with respect to
principal pivot transforms and pivotal method to its solution point of view.
For details see [3], [28], [4] and [29].

There so many methods are developed to solve a nonlinear complementarity
problem. see [31], [32],[12], [35]. Eaves and Saigal [7] formed an important
class of globally convergent methods for solving systems of non-linear
equations. Such methods have been used to constructively prove the existence
of solutions to many economic and engineering problems. The fundamental idea
of a homotopy continuation method is to solve a problem by tracing a certain
continuous path that leads to a solution of the problem. Thus, defining a
homotopy mapping that yields a finite continuation path plays an essential
role in a homotopy continuation method.

The paper is organized as follows. Section 2 presents some basic notations and
results. In section 3, we propose a new homotopy function to find the solution
of nonlinear complementarity problem. We construct a smooth and bounded
homotopy path to find the solution of the nonlinear complementarity problem as
the homotopy parameter $\lambda$ tends to $0$. To find the solution of
homotopy function we use modified homotopy continuation method to increase the
order of convergency of the algorithm. We also find the sign of the positive
tangent direction of the homotopy path. Finally, in section 4, we numerically
solve the oligopoly problem which is formulated by nonlinear complementarity
problem using the introduced homotopy function.

##  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. 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,\lambda)=$
$(1-\lambda)f(x)+\lambda g(x),$ $0\leq\lambda\leq 1.$ Then we consider the
homotopy equation $H(x,\lambda)=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,\lambda)=0$ varrying the values of $\lambda$ from $1$ to $0$. Kojima et
al showed that under some conditions nonlinear complementarity problem can be
solvable by homotopy continuation method. For details see [14], [13], [15],
[33].

Now we state some results which will be required in the next section.

######  Lemma 2.1:

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

######  Lemma 2.2:

(Parameterized Sard Theorem [34]) Let $V\subset R^{n},U\subset R^{m}$ be open
sets, and let $\phi:V\times U\to R^{k}$ be a $C^{\alpha}$ mapping, where
$\alpha>\text{max}\\{0,m-k\\}.$ If $0\in R^{k}$ is a regular value of $\phi,$
then for almost all $a\in V,0$ is a regular value of $\phi_{a}=\phi(a,.).$

######  Lemma 2.3:

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

######  Lemma 2.4:

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

######  Lemma 2.5:

[2] Consider the function $f:R^{n}\to R^{n}$ and the iterative method
$y^{k}=x^{k}-f^{\prime}(x^{k})^{-1}f(x^{k}),$
$z^{k}=x^{k}-2(f^{\prime}(y^{k})+f^{\prime}(x^{k}))^{-1}f(x^{k}),\
w^{k}=z^{k}-f^{\prime}(y^{k})^{-1}f(z^{k})$ has $5$th order convergence.

##  3 Main Results

Now we solve nonlinear complementarity problem by homotopy method. Consider
two positive numbers $m,n,$ such that $m$ is very large positive number and
$l$ is positive number, $l<<m.$ First we define
$\mathcal{R}_{(m)}=\\{(z,y,w_{1},w_{2},v_{1},v_{2})\in R_{++}^{n}\times
R_{++}^{n}\times R_{++}^{n}\times R_{++}^{n}\times R_{++}\times
R_{++}:m-(\overset{n}{\underset{i=1}{\sum}}(z+w_{1})_{i}+v_{2})>l,m-(\overset{n}{\underset{i=1}{\sum}}(y+w_{2})_{i}+v_{1})>l\\},\
$
$\mathcal{\bar{R}}_{(m)}=\\{(z,y,w_{1},w_{2},v_{1},v_{2})\in R_{+}^{n}\times
R_{+}^{n}\times R_{+}^{n}\times R_{+}^{n}\times R_{+}\times
R_{+}:m-(\overset{n}{\underset{i=1}{\sum}}(z+w_{1})_{i}+v_{2})\geq
l,m-(\overset{n}{\underset{i=1}{\sum}}(y+w_{2})_{i}+v_{1})\geq l\\}.$
Here $m$ is a predefined large number. We choose the initial point
$x^{(0)}=(z^{(0)},y^{(0)},{w_{1}}^{(0)},{w_{2}}^{(0)},{v_{1}}^{(0)},{v_{2}}^{(0)})\in\mathcal{R}_{(m)}$
such that
$A^{0}B^{0}-B^{0}{v_{2}}^{(0)}-{A^{(0)}{v_{1}}^{(0)}}=0,$
$l(B^{(0)}{v_{2}}^{(0)}-A^{(0)}{v_{1}}^{(0)})+lB^{(0)}(l-A^{(0)})+A^{(0)}{v_{1}}^{(0)}(A^{(0)}-{v_{2}}^{(0)})\neq
0,\\\
l(A^{(0)}{v_{1}}^{(0)}-B^{(0)}{v_{2}}^{(0)})+lA^{(0)}(l-B^{(0)})+B^{(0)}{v_{2}}^{(0)}(B^{(0)}-{v_{1}}^{(0)})\neq
0,\\\
l(B^{(0)}-l)\neq(A^{(0)}-{v_{2}}^{(0)}){v_{1}}^{(0)},l(A^{(0)}-l)\neq(B^{(0)}-{v_{1}}^{(0)}){v_{2}}^{(0)},$
where
$A=(m-(\overset{n}{\underset{i=1}{\sum}}(z+w_{1})_{i})),A^{(0)}=(m-(\overset{n}{\underset{i=1}{\sum}}(z^{(0)}+w_{1}^{(0)})_{i})),\\\
B=(m-(\overset{n}{\underset{i=1}{\sum}}(y+w_{2})_{i})),$
$B^{(0)}=m-(\overset{n}{\underset{i=1}{\sum}}(y^{(0)}+w_{2}^{(0)})_{i}).$
Now define the feasible region
$\mathcal{F}_{(m)}=\\{(z,y,w_{1},w_{2},v_{1},v_{2})\in{\mathcal{R}}_{(m)}:{v_{1}}\neq\frac{\lambda(A^{(0)}-{v_{2}}^{(0)}){v_{1}}^{(0)}}{l};\\\
{v_{2}}\neq\frac{\lambda(B^{(0)}-{v_{1}}^{(0)}){v_{2}}^{(0)}}{l}\ \forall\
\lambda\in(0,1)\\},$
$\mathcal{\bar{F}}_{(m)}=\\{(z,y,w_{1},w_{2},v_{1},v_{2})\in\mathcal{\bar{R}}_{(m)}:{v_{1}}\neq\frac{\lambda(A^{(0)}-{v_{2}}^{(0)}){v_{1}}^{(0)}}{l};{v_{2}}\neq\frac{\lambda(B^{(0)}-{v_{1}}^{(0)}){v_{2}}^{(0)}}{l}\
\forall\ \lambda\in(0,1)\\}.$
$\partial\mathcal{F}_{(m)}=\\{(z,y,w_{1},w_{2},v_{1},v_{2})\in\partial\mathcal{R}_{(m)}:{v_{1}}\neq\frac{\lambda(A^{(0)}-{v_{2}}^{(0)}){v_{1}}^{(0)}}{l};{v_{2}}\neq\frac{\lambda(B^{(0)}-{v_{1}}^{(0)}){v_{2}}^{(0)}}{l}\
\forall\ \lambda\in(0,1)\\},$ where $\partial\mathcal{R}_{(m)}$ is the
boundary of $\mathcal{\bar{R}}_{(m)}.$

Now we construct the homotopy function

$H(x,x^{(0)},\lambda)=\left[\begin{array}[]{c}(1-\lambda)(y-w_{1}+v_{1}e+J_{f}^{t}(z-w_{2}+v_{2}e))+\lambda(z-z^{(0)})\\\
W_{1}z-\lambda W_{1}^{(0)}z^{(0)}\\\ W_{2}y-\lambda W_{2}^{(0)}y^{(0)}\\\
y-(1-\lambda)f(z)-\lambda(y^{(0)})\\\
(m-\overset{n}{\underset{i=1}{\sum}}(z+w_{1})_{i}-v_{2})v_{1}-\lambda((m-\overset{n}{\underset{i=1}{\sum}}(z^{(0)}+w_{1}^{(0)})_{i}-v_{2}^{(0)})v_{1}^{(0)})\\\
(m-\overset{n}{\underset{i=1}{\sum}}(y+w_{2})_{i}-v_{1})v_{2}-\lambda((m-\overset{n}{\underset{i=1}{\sum}}(y^{(0)}+w_{2}^{(0)})_{i}-v_{1}^{(0)})v_{2}^{(0)})\\\
\end{array}\right]=0$ (3.2)

where $e=[1,1,\cdots,1]^{t},\ Z=\text{diag}(z);\ W_{1}=\text{diag}(w_{1});\
W_{2}=\text{diag}(w_{2});\ W_{1}^{(0)}=\text{diag}(w_{1}^{(0)});\
W_{2}^{(0)}=\text{diag}(w_{2}^{(0)});\
x=(z,y,w_{1},w_{2},v_{1},v_{2})\in\mathcal{\bar{F}}_{(m)};\
x^{(0)}=(z^{(0)},y^{(0)},{w_{1}}^{(0)},{w_{2}}^{(0)},{v_{1}}^{(0)},{v_{2}}^{(0)})\in\mathcal{F}_{(m)};\
\lambda\in(0,1]$ and $J_{f}$ is the jacobian of $f(z).$

######  Theorem 3.1:

For almost all initial points $x^{(0)}\in\mathcal{F}_{(m)},$ $0$ is a regular
value of the homotopy function $H:R^{4n+2}\times(0,1]\to R^{4n+2}$ and the
zero point set
$H_{x^{(0)}}^{-1}(0)=\\{(x,\lambda)\in\mathcal{F}_{(m)}\times(0,1]:H_{x^{(0)}}(x,\lambda)=0\\}$
contains a smooth curve $\Gamma_{x}^{(0)}$ starting from $(x^{(0)},1).$

###### Proof.

The jacobian matrix of the above homotopy function $H(x,x^{(0)},\lambda)$ is
denoted by $DH(x,x^{(0)},\lambda))$ and we have
$DH(x,x^{(0)},\lambda))=$$\left[\begin{array}[]{ccc}\frac{\partial{H(x,x^{(0)},\lambda)}}{\partial{x}}&\frac{\partial{H(x,x^{(0)},\lambda)}}{\partial{x^{(0)}}}&\frac{\partial{H(x,x^{(0)},\lambda)}}{\partial{\lambda}}\\\
\end{array}\right].$ For all $x^{(0)}\in\mathcal{F}_{1}$ and
$\lambda\in(0,1],$ we have
$\frac{\partial{H(x,x^{(0)},\lambda)}}{\partial{x^{(0)}}}=$$\begin{bmatrix}K_{1}&K_{2}\\\
K_{3}&K_{4}\\\ \end{bmatrix},$
where $K_{1}=\begin{bmatrix}-\lambda I&0&0&0\\\ -\lambda
W_{1}^{(0)}&0&-\lambda Z^{(0)}&0\\\ 0&-\lambda W_{2}^{(0)}&0&-\lambda
Y^{(0)}\\\ 0&-\lambda I&0&0\\\ \end{bmatrix},$ $K_{2}=\begin{bmatrix}0&0\\\
0&0\\\ 0&0\\\ 0&0\\\ \end{bmatrix},$

$K_{3}=\begin{bmatrix}\lambda v_{1}^{(0)}e^{t}&0&\lambda v_{1}^{(0)}e^{t}&0\\\
0&\lambda v_{2}^{(0)}e^{t}&0&\lambda v_{2}^{(0)}e^{t}\\\ \end{bmatrix},$

$K_{4}=\begin{bmatrix}-\lambda(m-\overset{n}{\underset{i=1}{\sum}}(z^{(0)}+w_{1}^{(0)})_{i}-v_{2}^{(0)})&\lambda
v_{1}^{(0)}\\\ \lambda
v_{2}^{(0)}&-\lambda(m-\overset{n}{\underset{i=1}{\sum}}(y^{(0)}+w_{2}^{(0)})_{i}-v_{1}^{(0)})\\\
\end{bmatrix},$
$Y^{(0)}=\text{diag}(y^{(0)}),Z^{(0)}=\text{diag}(z^{(0)})$,
$W_{1}^{(0)}=\text{diag}(w_{1}^{(0)})$,
$W_{2}^{(0)}=\text{diag}(w_{2}^{(0)})$.
$\text{det}(\frac{\partial{H}}{\partial{x^{(0)}}})=\lambda^{4n+2}((m-\overset{n}{\underset{i=1}{\sum}}(z^{(0)}+w_{1}^{(0)})_{i}-v_{2}^{(0)})(m-\overset{n}{\underset{i=1}{\sum}}(y^{(0)}+w_{2}^{(0)})_{i}-v_{1}^{(0)})-v_{1}^{(0)}v_{2}^{(0)})\prod_{i=1}^{n}z_{i}^{(0)}y_{i}^{(0)}\neq
0$ for $\lambda\in(0,1].$
Thus $DH(x,x^{(0)},\lambda))$ is of full row rank. Therefore, $0$ is a regular
value of $H(x,x^{(0)},\lambda)).$ By Lemmas 2.2 and 2.3, for almost all
$x^{(0)}\in\mathcal{F}_{(m)},$ $0$ is a regular value of
$H_{x^{(0)}}(x,\lambda)$ and $H_{x^{(0)}}^{-1}(0)$ consists of some smooth
curves and $H_{x^{(0)}}(x^{(0)},1)=0.$ Hence there must be a smooth curve
$\Gamma_{x}^{(0)}$ starting from $(x^{(0)},1).$

######  Theorem 3.2:

Let $\mathcal{F}_{(m)}$ be a nonempty set. For a given
$x^{(0)}\in\mathcal{F}_{(m)},$ if $0$ is a regular value of
$H(x,x^{(0)},\lambda),$ then $\Gamma_{x}^{(0)}$ is a bounded curve in
$\mathcal{\bar{F}}_{(m)}\times(0,1].$

###### Proof.

We have that $0$ is a regular value of $H(x,x^{(0)},\lambda)$ by theorem 3.1
and $\mathcal{F}_{(m)}$ be a nonempty set. It is clear that the set
$\mathcal{F}_{(m)}$ and $(0,1]$ is bounded. So there exists a sequence of
points
$\\{z^{k},y^{k},w_{1}^{k},w_{2}^{k},v_{1}^{k},v_{2}^{k},\lambda^{k}\\}\subset\Gamma_{x}^{(0)}\times(0,1],$
such that
$\lim\limits_{k\to\infty}z^{k}=\bar{z},\lim\limits_{k\to\infty}y^{k}=\bar{y},\lim\limits_{k\to\infty}w_{1}^{k}=\bar{w_{1}},\lim\limits_{k\to\infty}w_{2}^{k}=\bar{w_{2}},\lim\limits_{k\to\infty}v_{1}^{k}=\bar{v_{1}},\lim\limits_{k\to\infty}v_{2}^{k}=\bar{v_{2}},\lim\limits_{k\to\infty}\lambda^{k}=\bar{\lambda}.$
Hence $\Gamma_{x}^{(0)}$ is a bounded curve in
$\mathcal{\bar{F}}_{(m)}\times(0,1].$

Now we show the convergence of the homotopy function 3.2.

######  Theorem 3.3:

For
$x^{(0)}=(z^{(0)},y^{(0)},w_{1}^{(0)},w_{2}^{(0)},v_{1}^{(0)},v_{2}^{(0)})\in\mathcal{R}_{(m)}$
such that
$A^{0}B^{0}-B^{0}{v_{2}}^{(0)}-{A^{(0)}{v_{1}}^{(0)}}=0,$
$l(B^{(0)}{v_{2}}^{(0)}-A^{(0)}{v_{1}}^{(0)})+lB^{(0)}(l-A^{(0)})+A^{(0)}{v_{1}}^{(0)}(A^{(0)}-{v_{2}}^{(0)})\neq
0,\\\
l(A^{(0)}{v_{1}}^{(0)}-B^{(0)}{v_{2}}^{(0)})+lA^{(0)}(l-B^{(0)})+B^{(0)}{v_{2}}^{(0)}(B^{(0)}-{v_{1}}^{(0)})\neq
0,\\\ l(B^{(0)}-l)\neq(A^{(0)}-{v_{2}}^{(0)}){v_{1}}^{(0)},\\\
l(A^{(0)}-l)\neq(B^{(0)}-{v_{1}}^{(0)}){v_{2}}^{(0)},$
the homotopy equation finds a bounded smooth curve
$\Gamma_{x}^{(0)}\subset\mathcal{{F}}_{(m)}\times(0,1]$ which starts from
$(x^{(0)},1)$ and approaches the hyperplane at $\lambda\to 0.$ As $\lambda\to
0,$ the limit set
$\mathcal{L}\times\\{0\\}\subset\mathcal{\bar{F}}_{(m)}\times\\{0\\}$ of
$\Gamma_{x}^{(0)}$ is nonempty and every point in $\mathcal{L}$ is a solution
of the following system of equations:

$\displaystyle(y-w_{1}+v_{1}e+J_{f}^{t}(z-w_{2}+v_{2}e))=0$ (3.3)
$\displaystyle W_{1}z=0$ $\displaystyle W_{2}y=0$ $\displaystyle y-f(z)=0$
$\displaystyle(m-\overset{n}{\underset{i=1}{\sum}}(z+w_{1})_{i}-v_{2})v_{1}=0$
$\displaystyle(m-\overset{n}{\underset{i=1}{\sum}}(y+w_{2})_{i}-v_{1})v_{2}=0$

###### Proof.

Note that $\Gamma_{x}^{(0)}$ is diffeomorphic to a unit circle or a unit
interval $(0,1]$ in view of Lemma 2.4. As
$\frac{\partial{H(x,x^{(0)},1)}}{\partial{x^{(0)}}}$ is nonsingular,
$\Gamma_{x}^{(0)}$ is diffeomorphic to a unit interval $(0,1].$ Again
$\Gamma_{x}^{(0)}$ is a bounded smooth curve by the Theorem 3.2. Let
$(\bar{x},\bar{\lambda})$ be a limit point of $\Gamma_{x}^{(0)}.$ We consider
four cases:

  1. (i)

$(\bar{x},\bar{\lambda})\in\mathcal{{F}}_{(m)}\times\\{1\\}.$

  2. (ii)

$(\bar{x},\bar{\lambda})\in\partial{\mathcal{{F}}_{(m)}}\times\\{1\\}.$

  3. (iii)

$(\bar{x},\bar{\lambda})\in\partial{\mathcal{{F}}_{(m)}}\times(0,1).$

  4. (iv)

$(\bar{x},\bar{\lambda})\in\mathcal{\bar{F}}_{(m)}\times\\{0\\}.$

Suppose for case (i) the homotopy function 3.2 has solution $(\bar{x},1)$,
other than the initial solution $x^{(0)}$. As $\lambda\to 1$,
$\bar{y}=y^{(0)},\bar{z_{1}}=z_{1}^{0},\bar{z_{2}}=z_{2}^{0},\bar{v_{1}}\neq
0,\bar{v_{2}}\neq 0$. So for $\lambda\to 1$,
$(A-v_{2})\to(A^{(0)}-\bar{v_{2}}),(B-v_{1})\to(B^{(0)}-\bar{v_{1}})$. Hence
from homotopy function 3.2

$\displaystyle(A^{(0)}-\bar{v_{2}})\bar{v_{1}}=(A^{0}-{v_{2}}^{(0)}){v_{1}}^{(0)}$
(3.4)
$\displaystyle(B^{(0)}-\bar{v_{1}})\bar{v_{2}}=(B^{0}-{v_{1}}^{(0)}){v_{2}}^{(0)}$
(3.5)

From 3.5,
$\bar{v_{2}}=\frac{B^{0}-{v_{1}}^{(0)}}{B^{(0)}-\bar{v_{1}}}{v_{2}}^{(0)}.$
From equation 3.4
$(A^{(0)}-\frac{B^{0}-{v_{1}}^{(0)}}{B^{(0)}-\bar{v_{1}}}{v_{2}}^{(0)})\bar{v_{1}}=(A^{0}-{v_{2}}^{(0)}){v_{1}}^{(0)}.$
This implies that
$\bar{v_{1}}=\frac{-(B^{0}{v_{2}}^{(0)}-A^{0}{v_{1}}^{(0)}-A^{0}B^{0})\pm\sqrt{(A^{0}B^{0}-B^{0}{v_{2}}^{(0)}-A^{0}{v_{1}}^{(0)})^{2}}}{2A^{0}}$.
$\implies\bar{v_{1}}={v_{1}}^{(0)}\text{or}$
$\frac{A^{0}B^{0}-B^{0}{v_{2}}^{(0)}}{A^{(0)}}$. As
$\bar{v_{1}}=\frac{A^{0}B^{0}-B^{0}{v_{2}}^{(0)}}{A^{(0)}}={v_{1}}^{(0)}$ from
the condition of choosing the initial point $x^{(0)},$ the equation
$H_{x^{(0)}}(x,1)=0$ has only one solution $x^{(0)}\in\mathcal{{R}}_{(m)}.$
Hence the case $(i)$ is impossible.
In case $(ii)$ the homotopy equation 3.2 implies that
$\bar{y}=y^{(0)},\bar{z_{1}}=z_{1}^{0},\bar{z_{2}}=z_{2}^{0},\bar{v_{1}}\neq
0,\bar{v_{2}}\neq 0.$ So $(A-v_{2})\to(A^{(0)}-\bar{v_{2}})$ and
$(B-v_{1})\to(B^{(0)}-\bar{v_{1}})$ as $\lambda\to 1.$ From last two
components of homotopy equation 3.2 we have

$(A^{(0)}-\bar{v_{2}})\bar{v_{1}}=(A^{(0)}-{v_{2}}^{(0)}){v_{1}}^{(0)},(B^{(0)}-\bar{v_{1}})\bar{v_{2}}=(B^{(0)}-{v_{1}}^{(0)}){v_{2}}^{(0)}.$
(3.6)

Three cases may arise.
Case 1: Let $A^{(0)}-\bar{v_{2}}=l.$ From equation 3.6
$\bar{v_{1}}=\frac{(A^{(0)}-{v_{2}}^{(0)}){v_{1}}^{(0)}}{l},\
(B^{(0)}-\frac{(A^{(0)}-{v_{2}}^{(0)}){v_{1}}^{(0)}}{l})\bar{v_{2}}=(B^{(0)}-{v_{1}}^{(0)}){v_{2}}^{(0)}\\\
\implies\bar{v_{2}}=$
$\frac{l(B^{(0)}-{v_{1}}^{(0)}){v_{2}}^{(0)}}{lB^{(0)}-A^{(0)}{v_{1}}^{(0)}+{v_{1}}^{(0)}{v_{2}}^{(0)}}=A^{(0)}-l\\\
\implies
l(B^{(0)}{v_{2}}^{(0)}-A^{(0)}{v_{1}}^{(0)})+lB^{(0)}(l-A^{(0)})+A^{(0)}{v_{1}}^{(0)}(A^{(0)}-{v_{2}}^{(0)})=0,$
contradicts the choosing of initial point.
Case 2: Let $B^{(0)}-\bar{v_{1}}=l.$ From equation 3.6
$\bar{v_{2}}=\frac{(B^{(0)}-{v_{1}}^{(0)}){v_{2}}^{(0)}}{l},\
(A^{(0)}-\frac{(B^{(0)}-{v_{1}}^{(0)}){v_{2}}^{(0)}}{l})\bar{v_{1}}=(A^{(0)}-{v_{2}}^{(0)}){v_{1}}^{(0)}\\\
\implies\bar{v_{1}}=$
$\frac{l(A^{(0)}-{v_{2}}^{(0)}){v_{1}}^{(0)}}{lA^{(0)}-B^{(0)}{v_{2}}^{(0)}+{v_{1}}^{(0)}{v_{2}}^{(0)}}=B^{(0)}-l\\\
\implies
l(A^{(0)}{v_{1}}^{(0)}-B^{(0)}{v_{2}}^{(0)})+lA^{(0)}(l-B^{(0)})+B^{(0)}{v_{2}}^{(0)}(B^{(0)}-{v_{1}}^{(0)})=0,$
contradicts the choosing of initial point.
Case 3: Let $B^{(0)}-\bar{v_{1}}=l,\ A^{(0)}-\bar{v_{2}}=l.$ From equation 3.6
we have
$l\bar{v_{1}}=(A^{(0)}-{v_{2}}^{(0)}){v_{1}}^{(0)},l\bar{v_{2}}=(B^{(0)}-{v_{1}}^{(0)}){v_{2}}^{(0)}.\\\
\implies
l(B^{(0)}-l)=(A^{(0)}-{v_{2}}^{(0)}){v_{1}}^{(0)},l(A^{(0)}-l)=(B^{(0)}-{v_{1}}^{(0)}){v_{2}}^{(0)},$
contradicts the choosing of initial point.

In case $(iii)$ from homotopy equation 3.2 we have
$\bar{z}>0,\bar{y}>0,\bar{w_{1}}>0,\bar{w_{2}}>0.$ Three cases may arise.
Case1: Let $A-v_{2}\to\bar{A}-\bar{v_{2}}=l,$ where
$\bar{A}=(m-\overset{n}{\underset{i=1}{\sum}}(\bar{z}+\bar{w_{1}})_{i}.$ Then
from equation 3.6 we have
$\bar{v_{1}}=\frac{\bar{\lambda}(A^{(0)}-{v_{2}}^{(0)}){v_{1}}^{(0)}}{l},$
which contradicts that $\bar{v_{1}}\in\partial{\mathcal{{F}}_{(m)}}.$
Case2: Let $B-v_{1}\to\bar{B}-\bar{v_{1}}=l,$ where
$\bar{B}=(m-\overset{n}{\underset{i=1}{\sum}}(\bar{y}+\bar{w_{2}})_{i}.$ Then
from equation 3.6 we have
$\bar{v_{2}}=\frac{\bar{\lambda}(B^{(0)}-{v_{1}}^{(0)}){v_{2}}^{(0)}}{l},$
which contradicts that $\bar{v_{2}}\in\partial{\mathcal{{F}}_{(m)}}.$
Case3: Let $A-v_{2}\to\bar{A}-\bar{v_{2}}=l,B-v_{1}\to\bar{B}-\bar{v_{1}}=l.$
Then from equation 3.6 we have
$\bar{v_{1}}=\frac{\bar{\lambda}(A^{(0)}-{v_{2}}^{(0)}){v_{1}}^{(0)}}{l}$ and
$\bar{v_{2}}=\frac{\bar{\lambda}(B^{(0)}-{v_{1}}^{(0)}){v_{2}}^{(0)}}{l},$
which contradicts that $\bar{v_{1}}\in\partial{\mathcal{{F}}_{(m)}}$ and
$\bar{v_{2}}\in\partial{\mathcal{{F}}_{(m)}}.$
Therefore $(iv)$ is the only possible case. Hence
$\bar{x}=(\bar{z},\bar{y},\bar{w_{1}},\bar{w_{2}},\bar{v_{1}},\bar{v_{2}})$ is
a solution of the system of equations 3.3.

######  Remark 3.1:

From the homotopy function 3.2 as $\lambda\to 0$ we get
$\bar{y}-\bar{w}_{1}+\bar{J}_{f}^{t}(\bar{z}-\bar{w}_{2})=0,$
$\bar{y}=f(\bar{z})$ and $\bar{w}_{1i}\bar{z}_{i}=0,$
$\bar{w}_{2i}\bar{y}_{i}=0\ \forall i\in\\{1,2,\cdots n\\}$, where
$\bar{J}_{f}$ is the jacobian of $f(z)$ at the point $\bar{z}.$ Now
$\bar{w}_{1}$ and $\bar{w}_{2}$ can be decomposed as
$\bar{w}_{1}=\bar{y}-\Delta\bar{y}\geq 0$ and
$\bar{w}_{2}=\bar{z}-\Delta\bar{z}\geq 0.$ Now it is clear that
$\bar{y}_{i}\bar{z}_{i}=\Delta\bar{y}_{i}\bar{z}_{i}=\Delta\bar{z}_{i}\bar{y}_{i}\
\forall i\in\\{1,2,\cdots n\\}$ and
$\bar{J}_{f}^{t}\Delta\bar{z}+\Delta\bar{y}=0.$ This implies that
$(\bar{Z}\bar{J}_{f}^{t}+\bar{Y})\Delta\bar{z}=0,$ where
$\bar{Y}=$diag$(\bar{y})$ and $\bar{Z}=$diag$(\bar{z}).$

######  Theorem 3.4:

The component $\bar{z}$ of
$(\bar{z},\bar{y},\bar{w}_{1},\bar{w}_{2},0)\in\mathcal{L}\times\\{0\\}$ gives
the solution of the complementarity problem2.1 if and only if
$\Delta\bar{z}_{i}\Delta\bar{y}_{i}=0$ or $\bar{w}_{1i}+\bar{w}_{2i}>0\
\forall i\in\\{1,2,\cdots n\\}.$

###### Proof.

Suppose $\bar{z}\geq 0$ and $\bar{y}=f(\bar{z})\geq 0$ give the solution of
the complementarity problem2.1. Then $\bar{z}_{i}\bar{y}_{i}=0$ $\forall
i\in\\{1,2,\cdots n\\}.$ This implies that $\bar{z}_{i}=0$ or $\bar{y}_{i}=0$
$\forall i\in\\{1,2,\cdots n\\}.$ Now we consider the following cases.
Case 1: For atleast one $i\in\\{1,2,\cdots n\\},$ let
$\bar{z}_{i}>0,\bar{y}_{i}=0.$ In view of Remark 3.1, this implies that
$\Delta\bar{y}_{i}=0\implies\Delta\bar{z}_{i}\Delta\bar{y}_{i}=0.$
Case 2: For atleast one $i\in\\{1,2,\cdots n\\},$ let
$\bar{y}_{i}>0,\bar{z}_{i}=0.$ In view of 3.1, this implies that
$\Delta\bar{z}_{i}=0\implies\Delta\bar{z}_{i}\Delta\bar{y}_{i}=0.$
Case 3: For atleast one $i\in\\{1,2,\cdots n\\},$ let
$\bar{y}_{i}=0,\bar{z}_{i}=0.$ This implies that either
$\Delta\bar{y}_{i}\Delta\bar{z}_{i}=0$ or $\bar{w}_{1i}+\bar{w}_{2i}>0.$

Conversely, let consider $\Delta\bar{z}_{i}\Delta\bar{y}_{i}=0$ or
$\bar{w}_{1i}+\bar{w}_{2i}>0\ \forall i\in\\{1,2,\cdots n\\}.$ Let $\forall
i\in\\{1,2,\cdots n\\},\ \Delta\bar{z}_{i}\Delta\bar{y}_{i}=0$ implies either
$\Delta\bar{z}_{i}=0$ or $\Delta\bar{y}_{i}=0.$ This implies that
$\bar{y}_{i}\bar{z}_{i}=0\ \forall i\in\\{1,2,\cdots n\\}.$ Therefore
$\bar{y}$ and $\bar{z}$ give the solution of given complementarity problem
2.1. Let consider $\bar{w}_{1}+\bar{w}_{2}>0.$ Then three cases will arise.

Case 1: Let $\bar{w}_{1i}>0,\bar{w}_{2i}=0$ for atleast one $i\in\\{1,2,\cdots
n\\}.$ This implies that $\bar{z}_{i}=0$ and $\bar{y}_{i}\geq 0.$
Case 2: Let $\bar{w}_{1i}=0,\bar{w}_{2i}>0$ for atleast one $i\in\\{1,2,\cdots
n\\}.$ This implies that $\bar{z}_{i}\geq 0$ and $\bar{y}_{i}=0.$
Case 3: Let $\bar{w}_{1i}>0,\bar{w}_{2i}>0$ for atleast one $i\in\\{1,2,\cdots
n\\}.$ This implies that $\bar{z}_{i}=0$ and $\bar{y}_{i}=0.$
Considering the above three cases $\bar{z}$ and $\bar{y}$ solves the
compplementarity problem 2.1.

######  Theorem 3.5:

If the nonlinear function $f(z)$ is a $P_{0}$ function, then the component
$\bar{z}$ of
$(\bar{z},\bar{y},\bar{w}_{1},\bar{w}_{2},0)\in\mathcal{L}\times\\{0\\}$ gives
the solution of the nonlinear complementarity problem 2.1.

###### Proof.

Let $f(z)$ be a $P_{0}$ function. Then the jacobian matrix of the nonlinear
function at a point $z$, $\bar{J}_{f}$ is a $P_{0}$ matrix. Assume that the
component $\bar{z}$ of
$(\bar{z},\bar{y},\bar{w}_{1},\bar{w}_{2},0)\in\mathcal{L}\times\\{0\\}$ does
not give the solution of the nonlinear complementarity problem 2.1. Hence
$\Delta\bar{z}_{i}\Delta\bar{y}_{i}\neq 0$ and $\bar{w}_{1i}+\bar{w}_{2i}=0$
for atleast one $i$. Then $\Delta\bar{z}_{i}\neq 0,\Delta\bar{y}_{i}\neq
0,\bar{w}_{1i}=0,\bar{w}_{2i}=0.$ Now
$\bar{w}_{1i}=\bar{y}_{i}-\Delta\bar{y}_{i}=0$ and
$\Delta\bar{z}_{i}\Delta\bar{y}_{i}\neq 0$
$\implies\bar{y}_{i}=\Delta\bar{y}_{i}>0$. In similar way
$\bar{w}_{2i}=\bar{z}_{i}-\Delta\bar{z}_{i}=0$ and
$\Delta\bar{z}_{i}\Delta\bar{y}_{i}\neq 0$
$\implies\bar{z}_{i}=\Delta\bar{z}_{i}>0$. As
$(\bar{z},\bar{y},\bar{w}_{1},\bar{w}_{2},0)\in\mathcal{\bar{F}}_{(m)}\times\\{0\\}$,
$v_{1}=0,v_{2}=0$. From Equation 3.3,
$\Delta\bar{y}_{i}+({\bar{J}_{f}}^{t}\Delta\bar{z})_{i}=0.$ This implies that
$({\bar{J}_{f}}^{t}\Delta\bar{z})_{i}<0$ and also
$(\bar{z})_{i}({\bar{J}_{f}}^{t}\Delta\bar{z})_{i}<0.$ This contradicts that
$\bar{J}_{f}$ is a $P_{0}$-matrix. Therefore the component $\bar{z}$ of
$(\bar{z},\bar{y},\bar{w}_{1},\bar{w}_{2},0)\in\mathcal{L}\times\\{0\\}$ gives
the solution of the nonlinear complementarity problem 2.1.

######  Theorem 3.6:

Suppose the matrix $(\bar{Z}\bar{J}_{f}^{t}+\bar{Y})$ is nonsingular, where
$\bar{Y}=$diag$(\bar{y})$ and $\bar{Z}=$diag$(\bar{z}).$ Then $\bar{z}$ solves
the complementarity problem 2.1.

###### Proof.

Let $(\bar{Z}\bar{J}_{f}^{t}+\bar{Y})$ is nonsingular matrix. Now from remark
3.1 it is clear that $\Delta z=0.$ This implies that $\bar{z}$ solves the
complementarity problem 2.1.

######  Theorem 3.7:

If the jacobian matrix $\bar{J}_{f}$ has nonsingular principal minors, then
the component $\bar{z}$ of
$(\bar{z},\bar{y},\bar{w}_{1},\bar{w}_{2},0)\in\mathcal{L}\times\\{0\\}$
solves the nonlinear complementarity problem 2.1.

###### Proof.

Consider that the jacobian matrix $\bar{J}_{f}$ has nonsingular principal
minors. Then the transpose of the jacobian matrix $\bar{J}_{f}^{t}$ has
nonsingular principal minors. By theorem 3.6, if the matrix
$(\bar{Z}\bar{J}_{f}^{t}+\bar{Y})$ is nonsingular, then $\bar{z}$ solves the
nonlinear complementarity problem, where $\bar{Y}=\text{diag}(\bar{y}),$
$\bar{Z}=\text{diag}(\bar{z}).$ Let
$\tilde{\mathcal{A^{\prime}}}=\left[\begin{array}[]{cc}\bar{Y}&\bar{Z}\\\
-\bar{J}_{f}^{t}&I\\\ \end{array}\right]$. Then
$\det(\tilde{\mathcal{A^{\prime}}})=\det(\bar{Y}+\bar{Z}\bar{J}_{f}^{t})$.
Assume that the component $\bar{z}$ of
$(\bar{z},\bar{y},\bar{w}_{1},\bar{w}_{2},0)\in\mathcal{L}\times\\{0\\}$ is
not the solution of the nonlinear complementarity problem. Then there exists
atleast one $i$, such that $\bar{z}_{i}\bar{y}_{i}>0$. Without loss of
generality $\bar{y}$ and $\bar{z}$ can be represented as
$\bar{y}=\left[\begin{array}[]{c}\bar{y}_{p}\\\ \bar{y}_{q}\\\ \bar{o}_{r}\\\
\end{array}\right]$, $\bar{z}=\left[\begin{array}[]{c}\bar{o}_{p}\\\
\bar{z}_{q}\\\ \bar{z}_{r}\\\ \end{array}\right]$, where
$\bar{y}_{p}\in{R^{p}}_{++},\ \bar{y}_{q},\bar{x}_{q}\in{R^{q}}_{++},\
\bar{z}_{r}\in{R^{r}}_{++}$, $\bar{o}_{r}\in R^{r},\ \bar{o}_{p}\in R^{p}$ and
$\bar{o}_{r},\ \bar{o}_{p}$ are vectors with all zeros. Here
$(\bar{y}_{q})_{i}(\bar{z}_{q})_{i}>0$ and
$\bar{Y}=\text{diag}(\bar{y}),\bar{Z}=\text{diag}(\bar{z})$. Now we can
rewrite $\left[\begin{array}[]{cc}\bar{Y}&\bar{Z}\\\ -\bar{J}_{f}^{t}&I\\\
\end{array}\right]=\left[\begin{array}[]{cccccc}\bar{Y}_{p}&\bar{O}_{q}&\bar{O}_{r}&\bar{O}_{p}&\bar{O}_{q}&\bar{O}_{r}\\\
\bar{O}_{p}&\bar{Y}_{q}&\bar{O}_{r}&\bar{O}_{p}&\bar{Z}_{q}&\bar{O}_{r}\\\
\bar{O}_{p}&\bar{O}_{q}&\bar{O}_{r}&\bar{O}_{p}&\bar{O}_{q}&\bar{Z}_{r}\\\
M^{\prime}&B^{\prime}&C^{\prime}&I_{p}&\bar{O}_{q}&\bar{O}_{r}\\\
D^{\prime}&E^{\prime}&F^{\prime}&\bar{O}_{p}&I_{q}&\bar{O}_{r}\\\
G^{\prime}&H^{\prime}&K^{\prime}&\bar{O}_{p}&\bar{O}_{q}&I_{r}\\\
\end{array}\right]$, where
$-\bar{J}_{f}^{t}=\left[\begin{array}[]{ccc}M^{\prime}&B^{\prime}&C^{\prime}\\\
D^{\prime}&E^{\prime}&F^{\prime}\\\ G^{\prime}&H^{\prime}&K^{\prime}\\\
\end{array}\right]$,
$\bar{Y}=\left[\begin{array}[]{ccc}\bar{Y}_{p}&\bar{O}_{q}&\bar{O}_{r}\\\
\bar{O}_{p}&\bar{Y}_{q}&\bar{O}_{r}\\\ \bar{O}_{p}&\bar{O}_{q}&\bar{O}_{r}\\\
\end{array}\right]$,
$\bar{Z}=\left[\begin{array}[]{ccc}\bar{O}_{p}&\bar{O}_{q}&\bar{O}_{r}\\\
\bar{O}_{p}&\bar{Z}_{q}&\bar{O}_{r}\\\ \bar{O}_{p}&\bar{O}_{q}&\bar{Z}_{r}\\\
\end{array}\right]$, $\bar{Z}_{q}=\text{diag}(\bar{z}_{q})$,
$\bar{Z}_{r}=\text{diag}(\bar{z}_{r})$,
$\bar{Y}_{q}=\text{diag}(\bar{y}_{q})$,
$\bar{Y}_{p}=\text{diag}(\bar{y}_{p})$,
$\bar{O}_{p}=\text{diag}(\bar{o}_{p})$, $\bar{O}_{q}=\text{diag}(\bar{o}_{q})$
$\bar{O}_{r}=\text{diag}(\bar{o}_{r})$,
$M^{\prime},D^{\prime},G^{\prime},I_{p}\in R^{p\times p}$,
$B^{\prime},E^{\prime},H^{\prime},I_{q}\in R^{q\times q}$,
$C^{\prime},F^{\prime},K^{\prime},I_{r}\in R^{r\times r}$ and
$I_{p},I_{q},I_{r}$ are identity matrices. By elementary row operations we can
get
$\tilde{\mathcal{B^{\prime}}}=\left[\begin{array}[]{cccccc}I&\bar{O}_{q}&\bar{O}_{r}&\bar{O}_{p}&\bar{O}_{q}&\bar{O}_{r}\\\
\bar{O}_{p}&I&\bar{O}_{r}&\bar{O}_{p}&\bar{Z}_{q}{\bar{Y}_{q}}^{-1}&\bar{O}_{r}\\\
\bar{O}_{p}&\bar{O}_{q}&\bar{O}_{r}&\bar{O}_{p}&\bar{O}_{q}&I\\\
M^{\prime}&B^{\prime}&C^{\prime}&I&\bar{O}_{q}&\bar{O}_{r}\\\
D^{\prime}&E^{\prime}&F^{\prime}&\bar{O}_{p}&I&\bar{O}_{r}\\\
G^{\prime}&H^{\prime}&K^{\prime}&\bar{O}_{p}&\bar{O}_{q}&I\\\
\end{array}\right]$. By interchanging rows this matrix reduces to
$\tilde{\mathcal{C^{\prime}}}=$
$\left[\begin{array}[]{cccccc}I&\bar{O}_{q}&\bar{O}_{r}&\bar{O}_{p}&\bar{O}_{q}&\bar{O}_{r}\\\
\bar{O}_{p}&I&\bar{O}_{r}&\bar{O}_{p}&\bar{Z}_{q}{\bar{Y}_{q}}^{-1}&\bar{O}_{r}\\\
-G^{\prime}&-H^{\prime}&-K^{\prime}&\bar{O}_{p}&\bar{O}_{q}&\bar{O}_{r}\\\
M^{\prime}&B^{\prime}&C^{\prime}&I&\bar{O}_{q}&\bar{O}_{r}\\\
D^{\prime}&E^{\prime}&F^{\prime}&\bar{O}_{p}&I&\bar{O}_{r}\\\
G^{\prime}&H^{\prime}&K^{\prime}&\bar{O}_{p}&\bar{O}_{q}&I\\\
\end{array}\right]$.
Hence
$\det(\tilde{\mathcal{A^{\prime}}})=\det(\tilde{\mathcal{C^{\prime}}})=(-1)^{r}\det(K^{\prime})\neq
0$. Therefore by theorem 3.6 $\bar{z}$ solves the nonlinear complementarity
problem. This contradicts our assumption. Hence the component $\bar{z}$ of
$(\bar{z},\bar{y},\bar{w}_{1},\bar{w}_{2},0)\in\mathcal{L}\times\\{0\\}$ is
the solution of the nonlinear complementarity problem 2.1.

######  Remark 3.2:

Now we trace the homotopy path
$\Gamma_{x}^{(0)}\subset\mathcal{{F}}_{(m)}\times(0,1]$ from the initial point
$(x^{(0)},1)$ until $\lambda\to 0$ and find the solution of the given
complementarity problem 2.1 under some assumptions. Let $s$ denote the arc
length of $\Gamma_{x}^{(0)}$, we can parameterize the homotopy path
$\Gamma_{x}^{(0)}$ with respect to $s$ in the following form

$H_{x^{(0)}}(x(s),\lambda(s))=0,\ x(0)=x^{(0)},\ \lambda(0)=1.$ (3.7)

Now differentiating 3.7 with respect to $s$ we obtain the following system of
ordinary differential equations with given initial values[8]

$H^{\prime}_{x^{(0)}}(x(s),\lambda(s))\left[\begin{array}[]{c}\frac{dx}{ds}\\\
\frac{d\lambda}{ds}\\\ \end{array}\right]=0,\
\|(\frac{dx}{ds},\frac{d\lambda}{ds})\|=1,\ x(0)=x^{(0)},\ \lambda(0)=1,\
\frac{d\lambda}{ds}(0)<0,$ (3.8)

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

Now we use the homotopy continuation method with some modifications to trace
the homotopy path $\Gamma_{x}^{(0)}$ numerically. For details see ( cite our
homotopy paper)

###  3.1 Algorithm

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

$\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}.$

$\epsilon_{1}=10^{-9},\epsilon_{2}=10^{-6}.$ These are real numbers, used as
thresholds for $\lambda.$ If $\lambda$ achieves a value
$0\leq\lambda\leq\epsilon_{1},$ then the algorithm stops with an Acceptable
Solution. But, due to a stuck out for some specific reasons, if $\lambda$
achieves a value, such that, $\epsilon_{1}<\lambda\leq\epsilon_{2},$ the
algorithm stops, declaring that point as Probable Solution.
. Step 1: Set $\left[\begin{array}[]{c}x\\\ t\\\
\end{array}\right]=\left[\begin{array}[]{c}x^{(0)}\\\ 1\\\
\end{array}\right].$ Now calculate the constant $d^{(0)}=\det(\frac{\partial
H}{\partial x}(x^{(0)},\lambda_{0})).$ If $|d^{(0)}|\leq\epsilon,$ then stop
else go to step $2$.[$\epsilon\to 0,$ a threshold.]

Step 2: Set $c_{1}=c_{2}=0.$ Now calculate the constant $d=\det(\frac{\partial
H}{\partial x}(x,\lambda)).$ If $|d|\leq\epsilon,$ then stop else go to step
$3$.[$\epsilon\to 0,$ a threshold.]

Step 3: Determine the unit predictor direction $\tau^{(n)}$ by the following
method: If sign$(d)=-\text{sign}(d_{0}),$ then set $t_{d}=1-\lambda,$ else set
$t_{d}=-\lambda.$ Calculate $w_{d}=-t_{d}(\frac{\partial H}{\partial
x}(x,\lambda))^{-1}(\frac{\partial H}{\partial\lambda}(x,\lambda)),$
$\tau^{(n)}=$$\left[\begin{array}[]{c}x_{n}\\\ t_{n}\\\
\end{array}\right]=\frac{1}{\|x_{d},t_{d}\|}\left[\begin{array}[]{c}x_{d}\\\
t_{d}\\\ \end{array}\right],$ $\tau=\dfrac{|t_{d}|}{\|x_{d},t_{d}\|},$ where
$\|x_{d},t_{d}\|=\sqrt{x_{d}^{2}+t_{d}^{2}}.$ If $\tau\leq\eta_{1},$ then set
$c1=c1+1$ else reset $c1=0.$ If $c1<c0,$ then go to step $4$ else, if
$t_{d}\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(x)]^{t}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{x})\leq\mu(x),$ and
$x,\bar{x}\in\mathcal{F}_{(}m),$ where $\bar{x}$ is the solution of the
problem. $\mu(x)$ is taken as $[H_{0}(x)]^{t}[H_{0}(x)],$ where
$H_{0}(x)=\left[\begin{array}[]{c}(y-w_{1}+v_{1}e+J_{f}^{t}(z-w_{2}+v_{2}e))\\\
W_{1}z\\\ W_{2}y\\\ y-f(z)\\\
(m-\overset{n}{\underset{i=1}{\sum}}(z+w_{1})_{i}-v_{2})v_{1}\\\
(m-\overset{n}{\underset{i=1}{\sum}}(y+w_{2})_{i}-v_{1})v_{2}\\\
\end{array}\right].$

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

else if
$\gamma<0,\mu(x+\kappa_{1}^{k+1}x_{n})<\mu(x+\kappa_{1}^{k}x_{n}),x+\kappa_{1}^{k+1}x_{n}\in\mathcal{F}_{(}m),0<t+\kappa_{1}^{k+1}t_{n}<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
$t_{n}\leq\epsilon_{2},$ then stop with probable solution else stop.

Step 7: Compute the predictor and corrector point:
$\left[\begin{array}[]{c}x_{p}\\\ t_{p}\\\
\end{array}\right]=\left[\begin{array}[]{c}x\\\ t\\\
\end{array}\right]+\kappa_{1}^{k}\left[\begin{array}[]{c}x_{n}\\\ t_{n}\\\
\end{array}\right],$ $\left[\begin{array}[]{c}\bar{x}_{p}\\\ \bar{t}_{p}\\\
\end{array}\right]=\left[\begin{array}[]{c}x_{p}\\\ t_{p}\\\
\end{array}\right]-[J_{H}(x_{p},t_{p})^{+}H(x_{p},t_{p})],$ where
$[J_{H}(x_{p},t_{p})^{+}]$ is the Moore-Penrose Inverse. Now compute
$\left[\begin{array}[]{c}\tilde{x}_{p}\\\ \tilde{t}_{p}\\\
\end{array}\right]=$ $\left[\begin{array}[]{c}x_{p}\\\ t_{p}\\\
\end{array}\right]-2[(J_{H}(x_{p},t_{p})+J_{H}(\bar{x}_{p},\bar{t}_{p}))^{+}H(x_{p},t_{p})]$
and $r=\|H(x_{c},t_{c})\|.$ Then compute the next iteration
$\left[\begin{array}[]{c}x_{cc}\\\ t_{cc}\\\
\end{array}\right]=\left[\begin{array}[]{c}\bar{x}_{p}\\\ \bar{t}_{p}\\\
\end{array}\right]-2[J_{H}(\bar{x}_{p},\bar{t}_{p})+J_{H}(\tilde{x}_{p},\tilde{t}_{p})]^{+}H(\bar{x}_{p},\bar{t}_{p}).$
Then $\left[\begin{array}[]{c}x_{b}\\\ t_{b}\\\
\end{array}\right]=\left[\begin{array}[]{c}x_{cc}\\\ t_{cc}\\\
\end{array}\right]-J_{H}(x_{cc},t_{cc})^{+}H(x_{cc},t_{cc}).$ Repeating
$m_{0}$ times get the next iteration $\left[\begin{array}[]{c}x_{c}\\\
t_{c}\\\ \end{array}\right]=\left[\begin{array}[]{c}x_{b}\\\ t_{b}\\\
\end{array}\right].$ If $r\leq 1,0<t_{c}<1,$ and $x_{c}\in\mathcal{F}_{(}m),$
then jump to step $10$ else set $k=k-1$ and go to step $8$.

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

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

Step 10: Set $\left[\begin{array}[]{c}x\\\ t\\\ \end{array}\right]=$
$\left[\begin{array}[]{c}x_{c}\\\ t_{c}\\\ \end{array}\right].$ If
$t_{c}\leq\epsilon_{1},$ then stop with acceptable homotopy solution else set
$i=i+1$ and go to step $2.$

Note that $J_{H}(x,t)^{+}$ is the moore penrose inverse of the jacobian matrix
$J_{H}(x,t).$ That is
$J_{H}(x,t)^{+}=J_{H}(x,t)^{T}(J_{H}(x,t)J_{H}(x,t)^{T})^{-1}.$ The proposed
homotopy continuation method solves homotopy function by solving the initial
value problem with the following iterative process $J_{j}=$
$\left[\begin{array}[]{c}x_{p}\\\ t_{p}\\\
\end{array}\right]=\left[\begin{array}[]{c}x\\\ t\\\
\end{array}\right]+\kappa_{1}^{k}\left[\begin{array}[]{c}x_{n}\\\ t_{n}\\\
\end{array}\right],$
$T_{j}=[J_{H}(x_{p},t_{p})^{+}H(x_{p},t_{p})],$
$S_{j}=J_{j}-K_{j}=$ $\left[\begin{array}[]{c}\bar{x}_{p}\\\ \bar{t}_{p}\\\
\end{array}\right]$,
$TT_{j}=$ $\left[\begin{array}[]{c}x_{p}\\\ t_{p}\\\
\end{array}\right]-2[(J_{H}(x_{p},t_{p})+J_{H}(\bar{x}_{p},\bar{t}_{p}))^{+}H(x_{p},t_{p})]=\left[\begin{array}[]{c}\tilde{x}_{p}\\\
\tilde{t}_{p}\\\ \end{array}\right]$.
$SS_{j}=$ $\left[\begin{array}[]{c}x_{cc}\\\ t_{cc}\\\
\end{array}\right]=\left[\begin{array}[]{c}\bar{x}_{p}\\\ \bar{t}_{p}\\\
\end{array}\right]-2[J_{H}(\bar{x}_{p},\bar{t}_{p})+J_{H}(\tilde{x}_{p},\tilde{t}_{p})]^{+}H(\bar{x}_{p},\bar{t}_{p}),$
By this iterative process the proposed homotopy function achieves the order of
convergence as $7^{m}-1.$

######  Theorem 3.8:

Suppose that the homotopy function has derivative, which is lipschitz
continuous in a convex neighbourhood $\cal N$ of $c$ where $c$ is the solution
of the homotopy function $H(u,t)=0,$ whose Jacobian matrix is continuous and
nonsingular and bounded on $\cal N.$ Then the homotopy continuation method has
order $7^{m}-1.$

###### Proof.

By the Implicit Function Theorem ensures the existence of a unique continuous
solution $z(h)\in\cal N$ of $\dot{z}(h)=-\tilde{J}^{-1}\tilde{f},$ $z(0)=u$
and $h\in(-\delta,\delta),$ for some $\delta>0.$ Define
$\beta_{j}=\|z(h)-I_{j}(u,h)\|.$ From lemma 2.5 $\beta_{j}=O(h^{7^{j}}).$ Then
$\beta_{j+1}=\|z(h)-I_{j+1}\|\leq K{\beta_{j}}^{7}.$ Hence
$\beta_{j+1}=O(h^{7^{j+1}}).$ By induction method the modified homotopy
continuation method has convergency of order $7^{m}-1$

######  Theorem 3.9:

If the homotopy curve $\Gamma_{x}^{(0)}$ is smooth, then the positive tangent
direction $\tau^{(0)}$ at the initial point $x^{(0)}$ satisfies
sign($\det\left[\begin{array}[]{c}\frac{\partial H}{\partial
x\partial\lambda}(x^{(0)},1)\\\ \tau^{(0)^{t}}\\\ \end{array}\right]$)$<0.$

###### Proof.

From equation 3.2 we have $H(x,x^{(0)},\lambda)=\\\ $
$\left[\begin{array}[]{c}(1-\lambda)(y-w_{1}+v_{1}e+J_{f}^{t}(z-w_{2}+v_{2}e))+\lambda(z-z^{(0)})\\\
W_{1}z-\lambda W_{1}^{(0)}z^{(0)}\\\ W_{2}y-\lambda W_{2}^{(0)}y^{(0)}\\\
y-(1-\lambda)f(z)-\lambda(y^{(0)})\\\
(m-\overset{n}{\underset{i=1}{\sum}}(z+w_{1})_{i}-v_{2})v_{1}-\lambda((m-\overset{n}{\underset{i=1}{\sum}}(z^{(0)}+w_{1}^{(0)})_{i}-v_{2}^{(0)})v_{1}^{(0)})\\\
(m-\overset{n}{\underset{i=1}{\sum}}(y+w_{2})_{i}-v_{1})v_{2}-\lambda((m-\overset{n}{\underset{i=1}{\sum}}(y^{(0)}+w_{2}^{(0)})_{i}-v_{1}^{(0)})v_{2}^{(0)})\\\
\end{array}\right]=0.$
Now at the point $(x=x^{(0)},\lambda=1)$ the value of the partial derivative
is $\frac{\partial H}{\partial
x\partial\lambda}(x,\lambda)=\begin{bmatrix}K_{5}&K_{6}\end{bmatrix},$ where
$K_{5}=\begin{bmatrix}M^{\prime}&N^{\prime}\\\ \end{bmatrix},$
$M^{\prime}=\begin{bmatrix}I&0&0&0\\\ W_{1}^{(0)}&0&Z^{(0)}&0\\\
0&W_{2}^{(0)}&0&Y^{(0)}\\\ 0&I&0&0\\\
-v_{1}^{(0)}e^{t}&0&-v_{1}^{(0)}e^{t}&0\\\
0&-v_{2}^{(0)}e^{t}&0&-v_{2}^{(0)}e^{t}\\\ \end{bmatrix}$
and $N^{\prime}=\begin{bmatrix}0&0\\\ 0&0\\\ 0&0\\\ 0&0\\\
(m-\overset{n}{\underset{i=1}{\sum}}(z^{(0)}+w_{1}^{(0)})_{i}-v_{2}^{(0)})&-v_{1}^{(0)}\\\
-v_{2}^{(0)}&(m-\overset{n}{\underset{i=1}{\sum}}(y^{(0)}+w_{2}^{(0)})_{i}-v_{1}^{(0)})\\\
\end{bmatrix}.$
$K_{6}=\begin{bmatrix}A\\\ B\\\ C\\\ D\\\ E\\\ F\\\
\end{bmatrix},$$Y^{(0)}=\text{diag}(y^{(0)}),Z^{(0)}=\text{diag}(z^{(0)})$,
$W_{1}^{(0)}=\text{diag}(w_{1}^{(0)})$,
$W_{2}^{(0)}=\text{diag}(w_{2}^{(0)}),\\\
A=-[y^{(0)}-w_{1}^{(0)}+v_{1}^{(0)}e+J_{(f^{(0)})}^{t}(z^{(0)}-w_{2}^{(0)}+v_{2}^{(0)}e)],\\\
B=-W_{1}^{(0)}z^{(0)},C=-W_{2}^{(0)}y^{(0)},D=f(z^{(0)})-y^{(0)},E=-(m-\overset{n}{\underset{i=1}{\sum}}(z^{(0)}+w_{1}^{(0)})_{i}-v_{2}^{(0)})v_{1}^{(0)},F=-(m-\overset{n}{\underset{i=1}{\sum}}(y^{(0)}+w_{2}^{(0)})_{i}-v_{1}^{(0)})v_{2}^{(0)}.$
Let positive tangent direction be $\tau^{(0)}=\left[\begin{array}[]{c}t\\\
-1\end{array}\right]=\left[\begin{array}[]{c}(R^{(0)}_{1})^{(-1)}R_{2}^{(0)}\\\
-1\end{array}\right],$
where $R^{(0)}_{1}=\left[\begin{array}[]{cc}P^{\prime}&Q^{\prime}\\\
\end{array}\right],$ $P^{\prime}=\left[\begin{array}[]{cccc}I&0&0&0\\\
W_{1}^{(0)}&0&Z^{(0)}&0\\\ 0&W_{2}^{(0)}&0&Y^{(0)}\\\ 0&I&0&0\\\
-v_{1}^{(0)}e^{t}&0&-v_{1}^{(0)}e^{t}&0\\\
0&-v_{2}^{(0)}e^{t}&0&-v_{2}^{(0)}e^{t}\\\ \end{array}\right],$
$Q^{\prime}=\left[\begin{array}[]{cc}0&0\\\ 0&0\\\ 0&0\\\
(m-\overset{n}{\underset{i=1}{\sum}}(z^{(0)}+w_{1}^{(0)})_{i}-v_{2}^{(0)})&-v_{1}^{(0)}\\\
-v_{2}^{(0)}&(m-\overset{n}{\underset{i=1}{\sum}}(y^{(0)}+w_{2}^{(0)})_{i}-v_{1}^{(0)})\\\
\end{array}\right]$
and $R^{(0)}_{2}=\left[\begin{array}[]{c}A\\\ B\\\ C\\\ D\\\ E\\\ F\\\
\end{array}\right],$ where
$A=-[y^{(0)}-w_{1}^{(0)}+v_{1}^{(0)}e+J_{(f^{(0)})}^{t}(z^{(0)}-w_{2}^{(0)}+v_{2}^{(0)}e)],\\\
B=-W_{1}^{(0)}z^{(0)},C=-W_{2}^{(0)}y^{(0)},D=f(z^{(0)})-y^{(0)},\\\
E=-(m-\overset{n}{\underset{i=1}{\sum}}(z^{(0)}+w_{1}^{(0)})_{i}-v_{2}^{(0)})v_{1}^{(0)},F=-(m-\overset{n}{\underset{i=1}{\sum}}(y^{(0)}+w_{2}^{(0)})_{i}-v_{1}^{(0)})v_{2}^{(0)}.$
Here $\text{det}(R^{(0)}_{1})=$
$(\frac{(m-\overset{n}{\underset{i=1}{\sum}}(z^{(0)}+w_{1}^{(0)})_{i}-v_{2}^{(0)})(m-\overset{n}{\underset{i=1}{\sum}}(y^{(0)}+w_{2}^{(0)})_{i}-v_{1}^{(0)})}{v_{1}^{(0)}}-v_{2}^{(0)})(m-\overset{n}{\underset{i=1}{\sum}}(y^{(0)}+w_{2}^{(0)})_{i}-\\\
v_{1}^{(0)})\prod_{i=1}^{n}z_{i}^{(0)}y_{i}^{(0)}\neq 0.$ Therefore,
$\det\left[\begin{array}[]{c}\frac{\partial H}{\partial
x\partial\lambda}(y^{(0)},1)\\\ \tau^{(0)^{t}}\\\
\end{array}\right]$$=\det\left[\begin{array}[]{cc}R^{(0)}_{1}&R^{(0)}_{2}\\\
(R^{(0)}_{2})^{t}(R^{(0)}_{1})^{(-t)}&-1\\\ \end{array}\right]$
$=\det\left[\begin{array}[]{cc}R^{(0)}_{1}&R^{(0)}_{2}\\\
0&-1-(R^{(0)}_{2})^{t}(R^{(0)}_{1})^{(-t)}(R^{(0)}_{1})^{(-1)}R_{2}^{(0)}\\\
\end{array}\right]\\\
=\det(R^{(0)}_{1})\det(-1-(R^{(0)}_{2})^{t}(R^{(0)}_{1})^{(-t)}(R^{(0)}_{1})^{(-1)}R_{2}^{(0)})\\\
=-\det(R^{(0)}_{1})\det(1+(R^{(0)}_{2})^{t}(R^{(0)}_{1})^{(-t)}(R^{(0)}_{1})^{(-1)}R_{2}^{(0)})\\\
=-(\frac{(m-\overset{n}{\underset{i=1}{\sum}}(z^{(0)}+w_{1}^{(0)})_{i}-v_{2}^{(0)})(m-\overset{n}{\underset{i=1}{\sum}}(y^{(0)}+w_{2}^{(0)})_{i}-v_{1}^{(0)})}{v_{1}^{(0)}}-v_{2}^{(0)})(m-\overset{n}{\underset{i=1}{\sum}}(y^{(0)}+w_{2}^{(0)})_{i}-\\\
v_{1}^{(0)})\prod_{i=1}^{n}z_{i}^{(0)}y_{i}^{(0)}\det(1+(R^{(0)}_{2})^{t}(R^{(0)}_{1})^{(-t)}(R^{(0)}_{1})^{(-1)}R_{2}^{(0)})<0.$

##  4 Numerical Example

We consider the nonlinear complementarity form of the oligopolistic market
equilibrium problem and determine the equilibrium point with homotopy method.

Here we consider the oligopolistic market equilibrium problem operating under
the Nash equilibrium concept of noncooperative behaviour, as a problem of the
game theory.

To illustrate the use of the above three algorithms presented in the previous
section, a numerical example[9] is presented here. 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.9)

The demand curve is given by:

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

The parameters of the equation 4.9 for the five firms is 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

Now we solve this problem using the above algorithm.

To solve this problem using the homotopy method 3.2 with real valued parameter
$\lambda,$ we first take the initial point
$z^{(0)}=$$\left[\begin{array}[]{c}1\\\ 1\\\ 1\\\ 1\\\ 1\\\
\end{array}\right],$ $y^{(0)}=$$\left[\begin{array}[]{c}1\\\ 1\\\ 1\\\ 1\\\
1\\\ \end{array}\right],$ ${w_{1}}^{(0)}=$$\left[\begin{array}[]{c}1\\\ 1\\\
1\\\ 1\\\ 1\\\ \end{array}\right],$
${w_{2}}^{(0)}=$$\left[\begin{array}[]{c}1\\\ 1\\\ 1\\\ 1\\\ 1\\\
\end{array}\right],$ $v_{1}=0.001,$ $v_{2}=0.001,$ $\lambda_{0}=1.$ After $23$
iterations we get the result $\bar{z}=$$\left[\begin{array}[]{c}15.42931\\\
12.49858\\\ 9.663473\\\ 7.165094\\\ 5.132566\\\ \end{array}\right],$
$\bar{y}=$$\left[\begin{array}[]{c}0\\\ 0\\\ 0\\\ 0\\\ 0\\\
\end{array}\right],$ $\bar{w_{1}}=$$\left[\begin{array}[]{c}0\\\ 0\\\ 0\\\
0\\\ 0\\\ \end{array}\right],$
$\bar{w_{2}}=$$\left[\begin{array}[]{c}15.42931\\\ 12.49858\\\ 9.663473\\\
7.165094\\\ 5.132566\\\ \end{array}\right],$ $\bar{v_{1}}=0,$ $\bar{v_{2}}=0,$
$\bar{\lambda}=0.$

Figure 1: Solution Approach of Oligopoly Problem

##  5 Conclusion

In this study, we consider homotopy path to solve nonlinear complementarity
problem based on our newly introduced homotopy function by modified homotopy
continuation method. We find the positive tangent direction of the homotopy
path. We prove that the smooth curve for the proposed homotopy function is
bounded and convergent under some conditions related to initial points. An
oligopoly equilibrium problem is numerically solved by the proposed modified
homotopy continuation method. .

##  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] Alicia Cordero, José L Hueso, Eulalia Martínez, and Juan R Torregrosa.  Increasing the convergence order of an iterative method for nonlinear systems.  Applied Mathematics Letters, 25(12):2369–2374, 2012.
  * [3] AK Das.  Properties of some matrix classes based on principal pivot transform.  Annals of Operations Research, 243, 05 2014.
  * [4] 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.
  * [5] 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.
  * [6] 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.
  * [7] B Curtis Eaves and Romesh Saigal.  Homotopies for computation of fixed points on unbounded regions.  Mathematical Programming, 3(1):225–237, 1972.
  * [8] Xiaona Fan, Tingting Xu, Furong Gao, et al.  Solving nonlinear complementarity problem by a smoothing homotopy method.  Taiwanese Journal of Mathematics, 19(1):51–63, 2015.
  * [9] Patrick Harker.  A variational inequality approach for the determination of ohgopolistic market equilibrium.  Mathematical Programming, 30:105–111, 09 1984.
  * [10] R Jana, AK Das, and A. Dutta.  On hidden z-matrix and interior point algorithm.  OPSEARCH, 56, 09 2019.
  * [11] 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.
  * [12] Stepan Karamardian.  The nonlinear complementarity problem with applications, part 1.  Journal of Optimization Theory and Applications, 4(2):87–98, 1969\.
  * [13] Masakazu Kojima, Nimrod Megiddo, and Shinji Mizuno.  A general framework of continuation methods for complementarity problems.  Mathematics of Operations Research, 18(4):945–963, 1993.
  * [14] Masakazu Kojima, Nimrod Megiddo, and Toshihito Noma.  Homotopy continuation methods for nonlinear complementarity problems.  Mathematics of Operations Research, 16(4):754–774, 1991.
  * [15] Masakazu Kojima, Toshihito Noma, and Akiko Yoshise.  Global convergence in infeasible-interior-point algorithms.  Mathematical Programming, 65(1):43–72, 1994.
  * [16] 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.
  * [17] 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.
  * [18] 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.
  * [19] Prasenjit Mondal, S Sinha, S.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.
  * [20] S.K. Neogy, R Bapat, AK Das, and T Parthasarathy.  Mathematical Programming and Game Theory for Decision Making.  04 2008.
  * [21] S.K. Neogy, R. Bapat, AK Das, and T. Parthasarathy.  Mathematical programming and game theory for decision making.  11 2021.
  * [22] 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.
  * [23] 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.
  * [24] S.K. Neogy and AK Das.  On almost type classes of matrices with q-property.  Linear & Multilinear Algebra , 53:243–257, 07 2005.
  * [25] S.K. Neogy and AK Das.  Some properties of generalized positive subdefinite matrices.  SIAM J. Matrix Analysis Applications, 27:988–995, 01 2006.
  * [26] S.K. Neogy and AK Das.  On weak generalized positive subdefinite matrices and the linear complementarity problem.  Linear and Multilinear Algebra, 61, 07 2013.
  * [27] S.K. Neogy, AK Das, and R. Bapat.  Modeling, computation and optimization.  06 2022.
  * [28] 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.
  * [29] 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.
  * [30] S.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.
  * [31] Jong-Shi Pang.  Inexact newton methods for the nonlinear complementarity problem.  Mathematical Programming, 36(1):54–71, 1986.
  * [32] Jong-Shi Pang and Steven A Gabriel.  Ne/sqp: A robust algorithm for the nonlinear complementarity problem.  Mathematical Programming, 60(1-3):295–337, 1993.
  * [33] Paul Tseng.  An infeasible path-following method for monotone complementarity problems.  SIAM Journal on Optimization, 7(2):386–402, 1997.
  * [34] Xiuyu Wang and Xingwu Jiang.  A homotopy method for solving the horizontal linear complementarity problem.  Computational and Applied Mathematics, 33, 04 2013.
  * [35] Layne T Watson.  Solving the nonlinear complementarity problem by a homotopy method.  SIAM Journal on Control and Optimization, 17(1):36–46, 1979.
  * [36] 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.