Thanks very much, Shenghui. Very helpful. I am not sure I understand the issue regarding Set vs SetDelayed with regard to memory usage. It is clear that using Set means the assembly of the very large symbolic expression is done only in the definition, and is then used as the template for substitution each time the function is called. I have no doubt that is better with respect to speed. But I thought that in either case the memory used in each call should be released after the call is complete, and not accumulate. However, the method using LinearModelFit you give in your first post is vastly superior to what I was doing -- I will use that!

Kind regards, David

If you rewrite the equation of the a circle at {x0,y0} with radius R0 in the following form:

x^2 + y^2 = (x0^2 + y0^2 -R0^2)  + (2* x0) * x  +(2*y0) * y  (* z = c + a * x + b * y *)

where {x,y} are points on the circumference. You can use the linear model fit to find the circle

In[23]:= (reg={2 #1,2 #2,#2^2+#1^2}&@@@points)//Short
Out[23]//Short= {{124,1042,275285},{132,1030,269581},{132,1044,276840},<<50295>>,{1756,1016,1028948},{1760,968,1008656}}

In[28]:= lm=LinearModelFit[reg,{1,x,y},{x,y}]
Out[28]= FittedModel[-323255.+471.705 x+501.105 y]

Best model fit by least square:

In[25]:= bf=lm["BestFitParameters"]
Out[25]= {-323255.,471.705,501.105}
In[26]:= exp=(x-#2)^2+(y-#3)^2-#1-#2^2-#3^2&@@bf
Out[26]= -150357.+(-471.705+x)^2+(-501.105+y)^2

The code is modified from this stack exchange page