Approximate strictly positive solution to a linear set of equations?

You can use the new in M12 function QuadraticOptimization, which minimizes functions of the form:

$$\frac{1}{2} x . q. x + c . x$$

subject to linear constraints on $x$. So, the first step is to figure out what $q$ and $c$ are for your example. To do this we expand your expr, but before doing so, note that for a vector u and a matrix M we have:

$$u.M=M^T.u$$

Then, we can expand expr yielding:

$$(M.X-b).(M.X-b) = X.M^T.M.X-2 b.M.X+b.b$$

Hence we have:

q := 2 Transpose[M] . M
c := -2 b . M

Let's do a small example and compare to @Roman's answer. Setup:

nn = 10;
SeedRandom[10];
M = Table[RandomReal[{0, 100}], {i, 1, nn}, {j, 1, nn}];
b = Table[RandomReal[{0, 100}], {i, 1, nn}];
X = Array[x, nn];
expr = M . X - b;

Comparison:

r1 = X /. Last @ Minimize[{expr.expr, X > 0}, X]
r2 = QuadraticOptimization[{q, c}, {IdentityMatrix[nn], ConstantArray[0, nn]}]

{0.221992, 0.188374, 0.131969, 0., 0., 0., 0., 0., 0.0849646, 0.028771}

{0.221992, 0.188374, 0.131969, 0., 0., 0., 0., 0., 0.0849646, 0.028771}

where I used constraints that ensured that each vector element is nonnegative. For your larger example:

nn = 1000;
SeedRandom[1];
M = Table[RandomReal[{0, 100}], {i, 1, nn}, {j, 1, nn}];
b = Table[RandomReal[{0, 100}], {i, 1, nn}];

res = QuadraticOptimization[{q, c}, {IdentityMatrix[nn], ConstantArray[0, nn]}]; //AbsoluteTiming
res[[;;20]]

{1.48153, Null}

{0., 0., 0., 0.015694, 0., 0.00561439, 0., 0.0157487, 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0., 0.}

So, it finishes on the order of a couple seconds, as requested.

Perhaps you could use the non-negative least-squares algorithm (NNLS) of Lawson and Hanson, discussed here. See the links in the comments.

NNLS was ported to Mathematica by Michael Woodhams, and used here as NNLS[M,b]. I am running Mma 11.3.0 for 64-bit Linux.

nn = 1000;
SeedRandom[10];
M = Table[RandomReal[{0, 100}], {i, 1, nn}, {j, 1, nn}];
b = Table[RandomReal[{0, 100}], {i, 1, nn}];
ls = NNLS[M,b];

Timing was about 10 seconds for nn=1000.