a(a-b) = 0.5(a^2 - b^2 + (a-b)^2)
a(a+b) = 0.5(a^2 + b^2 + (a-b)^2)
a(3a-4b+2c) = 0.5(a^2 - b^2 + (2a-b)^2 - (2b-c)^2 + (a-2b+c)^2)

(3a - 4b + c)(9a - 24b + 22c - 8d + e)

polynomialToSquares[poly_, vars_] := Module[
  {basis = Plus[##]^2 & @@@ Tuples[vars, {2}] // Union, params}
  , params = Block[{x}, Array[Unique[x] &, Length@basis]]
  ; basis.params /. First@SolveAlways[poly == basis.params, vars]
 ]
polynomialToSquares[poly_] := polynomialToSquares[poly, Variables@poly]

poly = a (a + b);
polynomialToSquares[poly]
poly - % // Expand
(* (3 a^2)/2 + b^2/2 - 1/2 (a + b)^2 *)
(* 0 *)

poly = a (3 a - 4 b + 2 c);
polynomialToSquares[poly]
poly - % // Expand
(* 4 a^2 + 2 b^2 - 2 (a + b)^2 - c^2 + (a + c)^2 *)
(* 0 *)

SeedRandom[10];
poly = Union[Times @@@ Tuples[vars, {2}]].RandomInteger[{-20, 20}, 10]
polynomialToSquares[poly]
poly - % // Expand
(* 3 a^2 + 15 a b + 18 b^2 + 2 a c + 12 b c + c^2 + 10 a d - 9 b d + 20 c d - 13 d^2 *)
(* -((21 a^2)/2) + 9 b^2 + 15/2 (a + b)^2 - 16 c^2 + (a + c)^2 + 6 (b + c)^2 - (47 d^2)/2 + 5 (a + d)^2 - 9/2 (b + d)^2 + 10 (c + d)^2 *)
(* 0 *)

ClearAll@polynomialToSquares
polynomialToSquares[poly_, vars_] := Module[{
   monomialBasis = Times @@@ Tuples[vars, {2}] // Union
   , squaredBasis = Plus[##]^2 & @@@ Tuples[vars, {2}] // Union
   , tempVars = Block[{x}, Array[Unique[x] &, 2 Length@vars]]
   , basisTransformation
  }
  , basisTransformation = Normal@Last@CoefficientArrays[
      Expand@squaredBasis /. Thread[monomialBasis -> tempVars],  tempVars
     ]
  ; Last@CoefficientArrays[poly /. Thread[monomialBasis -> tempVars], tempVars].Inverse[basisTransformation].squaredBasis
 ]
polynomialToSquares[poly_] := polynomialToSquares[poly, Variables@poly]

basis = {a^2, b^2, (a + b)^2};

First@SolveAlways[a (a - b) == Array[$x, 3].basis, {a, b}]
(* {$x[1] -> 3/2, $x[3] -> -(1/2), $x[2] -> 1/2} *)

Array[$x, 3].basis - a (a - b) /. {$x[1] -> 3/2, $x[3] -> -(1/2), $x[2] -> 1/2} // Expand
(* 0 *)

basis = {a^2, b^2, c^2, (a + b)^2, (a + c)^2, (b + c)^2};

RandomSeed[10];
poly = RandomInteger[{-10, 10}, 9].(Times @@@ Tuples[{a, b, c}, {2}])
(* -2 a^2 - 2 a b - 7 b^2 - 2 a c - 16 b c - 3 c^2 *)

soln = First@SolveAlways[poly == Array[$x, 6].basis, {a, b, c}];
Array[$x, 6].basis - poly /. soln // Expand
(* 0 *)