In the community of remote sensing, nonlinear mixing models have recently received particular attention in hyperspectral image processing. In this paper, we present a novel nonlinear spectral unmixing method following the recent multilinear mixing model of [1], which includes an infinite number of terms related to interactions between different endmembers. The proposed unmixing method is unsupervised in the sense that the endmembers are estimated jointly with the abundances and other parameters of interest, i.e., the transition probability of undergoing further interactions. Non-negativity and sum-to one constraints are imposed on abundances while only nonnegativity is considered for endmembers. The resulting unmixing problem is formulated as a constrained nonlinear optimization problem, which is solved by a block coordinate descent strategy, consisting of updating the endmembers, abundances and transition probability iteratively. The proposed method is evaluated and compared with linear unmixing methods for synthetic and real hyperspectral datasets acquired by the AVIRIS sensor. The advantage of using non-linear unmixing as opposed to linear unmixing is clearly shown in these examples.