We propose a new robust variation of the tensor decomposition known as the multi-linear generalized singular value decomposition (ML-GSVD), and demonstrate its effectiveness in the design of joint transmit (TX) and receive (RX) beamforming (BF) for both the downlink (DL) and the uplink (UL) of cell-free massive multiple-input multiple-output (CF-mMIMO) systems. The new technique, dubbed the robust orthogonality-enforcing ML-GSVD (ROEML-GSVD), mitigates the degradation in the DL and UL BF design due to the tensor decomposition error and the channel state information (CSI) imperfection by exploiting their independence. Numerical results demonstrate that the proposed method is highly effective, leading not only to a significant reduction on the tensor decomposition error compared to state-of-the-art (SotA) methods, but also to the improvement in the data rate achieved by the subsequent beamformers in heterogeneous scenarios where users have distinct numbers of antennas.INDEX TERMS Tensor decomposition, transmit and receive beamforming, heterogeneous cell-free massive MIMO systems.