You have to make more precise what what you mean by “approximately equal”. Typically, you would minimize \Vert x_1 A - x_2 B - x_3 C - D \Vert in some norm \Vert \cdots \Vert, but which norm do yo want?
@albheim’s solution suggested the L_2 (Euclidean) “vectorized” (Frobenius) norm, which makes this a constrained least-squares problem. (If you let x_3 = 1 - x_1 - x_2 to eliminate the equality constraint, then it is bound-constrained least-squares; see also this thread: Suggestions needed for bound-constrained least squares solver). If you use the L_1 or L_\infty vectorized norms, then it is equivalent to an LP (linear programming) problem. In any of these cases it is convex, and so it can be solved by lots of algorithms (e.g. via Convex.jl, JuMP.jl, StructuredOptimization.jl, …). With only 3 unknowns (or 2, if you eliminate x_3) and 10^6 equations, it’s not a particularly big problem, so lots of algorithms should be fine.