Actually you need a constraint for every data point. I can illustrate the approach with data that can be fitted using Fit. First generate the data and fit it to a straight line.
In[3]:= data = Table[{x, 2*x + 3 + RandomReal[{-.25, .25}]}, {x, 5}]
Out[3]= {{1, 4.86828}, {2, 7.17623}, {3, 9.03998}, {4, 11.2155}, {5,
13.1595}}
In[4]:= Fit[data, {1, x}, x]
Out[4]= 2.90539 + 2.06216 x
Now do the fit by minimizing the sum of the squares of the error terms.
In[5]:= {xvals, yvals} = Transpose[data]
Out[5]= {{1, 2, 3, 4, 5}, {4.86828, 7.17623, 9.03998, 11.2155,
13.1595}}
In[8]:= errterms = yvals - m*xvals - b
Out[8]= {4.86828 - b - m, 7.17623 - b - 2 m, 9.03998 - b - 3 m,
11.2155 - b - 4 m, 13.1595 - b - 5 m}
In[9]:= FindMinimum[errterms.errterms, {{m, 2}, {b, 3}}]
Out[9]= {0.0410029, {m -> 2.06216, b -> 2.90539}}
Now do it by defining constraints from the fitting equation. It is not implicit but the method would work if it were.
In[10]:= errtermsc = Table[yfit[i] - yvals[[i]], {i, Length[data]}]
Out[10]= {-4.86828 + yfit[1], -7.17623 + yfit[2], -9.03998 +
yfit[3], -11.2155 + yfit[4], -13.1595 + yfit[5]}
In[11]:= cons = Table[yfit[i] == m*xvals[[i]] + b, {i, Length[data]}]
Out[11]= {yfit[1] == b + m, yfit[2] == b + 2 m, yfit[3] == b + 3 m,
yfit[4] == b + 4 m, yfit[5] == b + 5 m}
In[12]:= vars =
Join[{{m, 2}, {b, 3}},
Table[{yfit[i], yvals[[i]]}, {i, Length[data]}]]
Out[12]= {{m, 2}, {b, 3}, {yfit[1], 4.86828}, {yfit[2],
7.17623}, {yfit[3], 9.03998}, {yfit[4], 11.2155}, {yfit[5], 13.1595}}
In[13]:= FindMinimum[{errtermsc.errtermsc, cons}, vars]
Out[13]= {0.0410029, {m -> 2.06216, b -> 2.90539, yfit[1] -> 4.96756,
yfit[2] -> 7.02972, yfit[3] -> 9.09189, yfit[4] -> 11.1541,
yfit[5] -> 13.2162}}
