The process of assembling isogeometric Galerkin matrices arising from hierarchical B-spline (HB-spline) discretizations is a topic of active research in view of the associated computational issues, especially when the dimension and the polynomial degree of the basis increase. To address this challenge, different approaches have been investigated. These include specialized quadrature rules, which reduce the number of evaluation points necessary for integration, and efficient Bézier extraction operators, which suitably combine isogeometric analysis with finite element codes. In addition, new assembly methods that go beyond classical element-wise algorithms leveraging the structure of the underlying basis were also proposed. In this talk, we propose a novel representation of HB-spline system matrices as block-wise Hadamard products, obtained through univariate integrals. The Hadamard format may be seen as a generalization of the Kronecker structure present in the tensor-product case. We use dedicated data structures to manipulate the sparse tensors involved in the assembly process, in order to reduce the memory footprint and the computational complexity of the method. A selection of numerical experiments will illustrate the effectiveness of our approach, implemented in the C++ library G+Smo.