In the present work, numerical methods for the solution of multi-linear system are presented. Most large-scale multi-linear solvers rely on either the alternating leastsquares or low rank Krylov methods. The approach we use to develop our methods lies somehow in between and can be considered as a generalisation of an alternated direction method. Given the multi-linear operator in the form of a sum of Kronecker product of matrices, we solve at each iteration a linear system for each summand. The approximate solution is then defined to be the best linear combination of these solutions, as well as the previous solution and the residual. Some convergence results are proved. Numerical experiments on two problems arising from parametric PDEs show the effectiveness of the proposed method
