Content Based Image Retrieval (CBIR) systems that search similar images in a large database are attracting more and more research interests recently, and have been applied to medical image characterization for expert's experience sharing. One challenging task in CBIR is how to extract features for effective image representation. Therein sparse coding technique has been proven to be an effective way to learn inherent structure features for image analysis. However, it is necessary to first vectorize the 2- or 3-dimensional spatial structure for analysis with sparse coding, and then destroy the spatial relation of nearby voxels. In this study, we propose a multilinear sparse coding method to learn features from multi-dimensional medical images. We regard high dimensional local structures as tensors and propose a K-CP (CANDECOMP/PARAFAC) algorithm to learn a tensor dictionary in an iterative way. With the learned tensor dictionary, sparse coefficients of tensor local structures are calculated by multilinear orthogonal matching pursuit (MOMP) algorithm, which is an extended multilinear version of the conventional linear OMP. The proposed multilinear sparse coding method is prospected to be more efficient and effective for inherent feature extraction compared with conventional linear methods. The proposed method is applied to a CBIR system for retrieval of focal liver lesions (FLLs) using a medical database consisting of contrast-enhanced multi-phase computer-tomography (CT) images. Experiments show that the constructed CBIR with multilinear sparse coding method can achieve promising retrieval performance.