We propose a novel numerical method for fast and accurate evaluation of the exchange part of the Fock operator in the Hartree-Fock equation which is a (nonlocal) integral operator. Usually, this challenging computational problem is solved by analytical evaluation of two-electron integrals using the “analytically separable” Galerkin basis functions, like Gaussians. Instead, we employ the agglomerated “grey-box” numerical computation of the corresponding six-dimensional integrals in the tensor-structured format which does not require analytical separability of the basis set. The point of our method is a low-rank tensor representation of arising functions and operators on an n×n×n Cartesian grid and the implementation of the corresponding multi-linear algebraic operations in the tensor product format. Linear scaling of the tensor operations, including the 3D convolution product, with respect to the one-dimension grid size n enables computations on huge 3D Cartesian grids thus providing the required high accuracy. The presented algorithm for evaluation of the exchange operator and a recent tensor method for the computation of the Coulomb matrix are the main building blocks in the numerical solution of the Hartree-Fock equation by the tensor-structured methods. These methods provide a new tool for algebraic optimization of the Galerkin basis in the case of large molecules.