###### Abstract
This formulation of the problem excludes the possibility of using more specialized linear algebra methods, which are based on a special structure of the stiffness matrix. In computational mathematics, numerical methods for evolution problems are naturally based on the best available solvers for the corresponding elliptic boundary value problems. In the case of standard (local) diffusion, the stiffness matrix is sparse, and the developed near-optimal implicit methods for parabolic problems use fast preconditioned iterative solution methods (such as algebraic multigrid) for the emergent linear systems. The situation changes substantially when fractional diffusion is involved. The matrices are dense, and the concept of iterative solvers is not applicable in the general case. Furthermore, it turns out that there are specific difficulties depending on the definition of the fractional Laplacian used. The history of the development of efficient (near-optimal) numerical methods for space-fractional diffusion problems in multidimensional domains of general geometry is less than ten years. And despite the considerable number of publications in this field, there are few results of this type for parabolic problems. Without pretending to be exhaustive, we will briefly touch on some of them. The case of spectral fractional diffusion is considered in [10; 11; 12; 13; 14]. The common feature of these results is that they can all be interpreted in terms of certain rational approximations. A further step in the development of this approach, significantly expanding the scope of the investigated problems, can be found in [15; 16; 17; 18]. These works mainly focus on convergence analysis, with less discussion of algorithmic implementation and computational complexity. The same applies to [19], where finite element approximations for integral-fractional evolution problems are analysed. This paper is thematically closest to the parabolic equation we are considering. Here, we use a rather different pure algebraic approach. Hierarchical semi-separable (HSS) compression is applied to approximately solve the arising linear systems. Earlier results based on HHS compression have been published in [20], where the case of a uniform time step is experimentally analysed. In this paper, we implement an adaptive forward-backward Euler scheme for discretization in time. The aim is to automatically adapt the method to the evolution of the regularity of the solution. The obtained composite algorithm is analysed. The block representation of HSS compression is used to upgrade the HSS solver. This allows us to efficiently handle different diagonally perturbed transition matrices corresponding to changing time steps. The main contributions of the paper are as follows. The method is fully constructive, leading to a well-structured algorithmic implementation. The theoretical estimate of the computational complexity clearly shows the efficiency of the adaptive time stepping in combination with the HSS solver. The presented numerical experiments support the theoretical results. Both sequential and parallel scalability and efficiency are analyzed. The conclusion is that a convincing proof of concept of the proposed method and algorithm is presented. In short, the fundamental novelties of the results that will be presented are the complex solutions obtained as a result of the developed composite algorithm. The time discretization is optimized, where the variable time step size adapts to the dynamics in the behavior of the solution. The hidden structure of the dense fractional diffusion stiffness matrix is shown to fit well with the stochastic HSS compression approach, thereby allowing attractive parallel scalability to be reached. The rest of the paper is organized as follows. The fractional-in-space parabolic problem is presented in the next section, which concludes with the adaptive forward-backward Euler scheme for time discretization. Section 3 is devoted to hierarchical semi-separable compression.