  This paper addresses multilinear systems of equations which arise in various
applications such as data mining and numerical partial differential equations.
When the multilinear system under consideration involves a nonsingular
$\mathcal{M}$-tensor and a nonnegative right-hand side vector, it may have
multiple nonnegative solutions. In this paper, we propose an algorithm which
can always preserve the nonnegativity of solutions. Theoretically, we show that
the sequence generated by the proposed algorithm is a nonnegative decreasing
sequence and converges to a nonnegative solution of the system. Numerical
results further support the novelty of the proposed method. Particularly, when
some elements of the right-hand side vector are zeros, the proposed algorithm
works well while existing state-of-the-art solvers may not produce a
nonnegative solution.
