This paper describes mathematical procedure for designing hexagonal systolic arrays that implement fault-tolerant matrix multiplication. Fault-tolerance is achieved by introducing redundancy at algorithm level by defining three equivalent algorithms with disjoint index spaces. The essence of the proposed method is based on mapping data dependency graph that corresponds to the matrix multiplication algorithm, by an appropriate epimorphism, into a graph with desired properties. Since there is a 1:1 correspondence between the algorithm and it's graph representation, all transformations performed on the graph directly affect the algorithm. Chosen epimorphism depends on the projection direction vector $\overrightarrow{\mu}=[\mu_1 \; \mu_2 \; \mu_3]^T$ and enables obtaining hexagonal arrays with optimal number of processing elements (PEs) for the given matrix dimensions, which realizes fault-tolerant matrix multiplication for the shortest possible time for that number of PEs. The proposed procedure is formally described by explicit formulas and can be used as a software tool for automatic synthesis of fault-tolerant arrays.