Linear algebraic expressions are the essence of many computationally intensive problems, including scientific simulations and machine learning applications. However, translating high-level formulations of these expressions to efficient machine-level representations is far from trivial: developers should be assisted by automatic optimization tools so that they can focus their attention on high-level problems, rather than low-level details. The tractability of these optimizations is highly dependent on the choice of the primitive constructs in terms of which the computations are to be expressed. In this work we propose to describe operations on multi-dimensional arrays using a selection of higher-order functions, inspired by functional programming, and we present rewrite rules for these such that they can be automatically optimized for modern hierarchical and heterogeneous architectures. Using this formalism we systematically construct and analyse different subdivisions and permutations of the dense matrix multiplication problem.