SummaryAn algorithm to generate samples with approximate first‐order, second‐order, third‐order, and fourth‐order moments is presented by extending the Cholesky matrix decomposition to a Cholesky tensor decomposition of an arbitrary order. The tensor decomposition of the first‐order, second‐order, third‐order, and fourth‐order objective moments generates a non‐linear system of equations. The algorithm solves these equations by numerical methods. The results show that the optimization algorithm delivers samples with an approximate residual error of less than 1016 between the components of the objective and the sample moments. The algorithm is extended for a n‐th‐order approximate tensor moment version, and simulations of non‐normal samples replicated from distributions with asymmetries and heavy tails are presented. An application for sensitivity analysis of portfolio risk assessment with Value‐at‐Risk (VaR) is provided. A comparison with previous methods available in the literature suggests that the methodology proposed reduces the error of the objective moments in the generated samples. Copyright © 2016 John Wiley & Sons, Ltd.