Supervised classification and representation learning are two widely usedclasses of methods to analyze multivariate images. Although complementary,these methods have been scarcely considered jointly in a hierarchical modeling.In this paper, a method coupling these two approaches is designed using amatrix cofactorization formulation. Each task is modeled as a factorizationmatrix problem and a term relating both coding matrices is then introduced todrive an appropriate coupling. The link can be interpreted as a clusteringoperation over a low-dimensional representation vectors. The attributionvectors of the clustering are then used as features vectors for theclassification task, i.e., the coding vectors of the correspondingfactorization problem. A proximal gradient descent algorithm, ensuringconvergence to a critical point of the objective function, is then derived tosolve the resulting non-convex non-smooth optimization problem. An evaluationof the proposed method is finally conducted both on synthetic and real data inthe specific context of hyperspectral image interpretation, unifying twostandard analysis techniques, namely unmixing and classification.