Can I use symbolic tensors for simple linear algebra and calculus?

Question

As an example, I'd like to calculate this symbolic derivative:

$\frac{\partial}{\partial b}(x-A.b)^{\mathsf{T}}.(x-A.b)$

where $x$ and $b$ are vectors and $A$ is a matrix.

What I tried is this:

Clear[x, a, b, n, m]    
$Assumptions = {x ∈ Vectors[n, Reals], a ∈ Matrices[{n, m}], b ∈ Vectors[m, Reals]}

expr = TensorReduce[(x - a.b).(x - a.b)]

The value returned by TensorReduce looks ok, in the sense that TensorRank and TensorDimensions return 0 and {} as expected for a scalar, and when I use replacement rules to replace x, a, b with proper-sized arrays, I get the expected result.

However, when I write:

D[expr, b]

I get the expression: -2 x.a.1 + a.1.a.b + a.b.a.1. Looking at this expression, I can sort of guess what the derivative should look like, but I can't use it in code like that. For example, I can't use TensorDimensions to get the derivative's dimension, and I can't use replacement rules to insert arrays in place of the symbolic tensors.

ADD:

I've also tried using TensorProduct/TensorContract instead of Dot:

xmab = TensorReduce[x - TensorContract[TensorProduct[a, b], {{2, 3}}]];    
expr = TensorReduce[
   TensorContract[TensorProduct[xmab, xmab], {{1, 2}}]]; 

This gives a sensible looking expr, but I think D doesn't work with TensorProduct/TensorContract:

D[expr, {b, 1}]

-2 (TensorContract^(1,{{0,0},{0,0}}))[a[TensorProduct]b[TensorProduct]x,{{1,4},{2,3}}] (TensorProduct^(0,1,0))[a,b,x]+(TensorContract^(1,{{0,0},{0,0},{0,0}}))[a[TensorProduct]a[TensorProduct]b[TensorProduct]b,{{1,3},{2,6},{4,5}}] ((TensorProduct^(0,0,0,1))[a,a,b,b]+(TensorProduct^(0,0,1,0))[a,a,b,b])

Answer 1

Finally, in version 14.1, VectorSymbol and MatrixSymbol are introduced:

x = VectorSymbol["x", n];
A = MatrixSymbol["A", {n, m}];
b = VectorSymbol["b", m];

rst = D[(x - A . b) . (x - A . b), b]

Let's test if the result is right:

bval = {b1, b2, b3};
xval = RandomInteger[{-10, 10}, {2}];
Aval = RandomInteger[{-10, 10}, {2, 3}];

Block[{b = bval, x = xval, A = Aval}, D[(x - A . b) . (x - A . b), {b}]] == 
   rst /. {b -> bval, x -> xval, A -> Aval} // Simplify
(* True *)