最小二乗法によって点群から直線のフィッテイングができるように，同じような手法で円のフィッティングもできないでしょうか？探してみると以下の素晴らしい記事がありました．
一般式による最小二乗法（円の最小二乗法）
記事によると，円の中心座標$(a,b)$，円の半径$r$としたときの円の一般式
(x-a)^{2}+(y-b)^{2}=r^{2}
の両辺の差分をとり，その二乗の総和
\sum^{n}_{i=1}((x_{i}-a)^{2}+(y_{i}-b)^{2}-r^{2})^{2} ⇒ min
が最小となる$a,b,r$を定めることによって円フィッテングが可能となる．というものです．
これ以上コピペしても仕方ありませんので，以降の解き方に関しましては是非リンク先を参照ください．
漠然と勾配降下法等で比較的大きな計算量が必要になると予想していたのですが，どうやら解析的に解くことができ，非常に軽量のようです．
すでに実装し公開して下さってる方もいらっしゃいました．コピペですぐ動作確認できますので，ぜひ試してみてください．
最小二乗法による点群の円フィッティングMATLAB＆Pythonサンプルプログラム
最小二乗法？
ところで，最初のリンク先のコメント欄で興味深い議論がありました．
この方法ですと、たとえば、点(3,0),(5,0),(0,3),(0,5),(-3,0),(-5,0),(0,-3),(0,-5)をフィッティングさせると、中心(0,0)、半径sqrt(17)の円になります。
本来は、中心(0,0)、半径4になることを期待していると思うのですが違うのでしょうか。
図示してみます．
図に示した点群からのフィッティングを行えば，感覚的にはそれぞれ距離の近い二点間の中心を円周が通ることになり，即ち半径4$=\sqrt{16}$の円となりそうですが，計算結果は$\sqrt{17}$となるようです．
手計算で確認してみます．距離が2となる同一軸上の二点$(3,0),(5,0)$で誤差の総和を計算してみると
(3^{2}-\sqrt{16}^{2})^{2}+(5^{2}-\sqrt{16}^{2})^{2}=130
(3^{2}-\sqrt{17}^{2})^{2}+(5^{2}-\sqrt{17}^{2})^{2}=128
確かに$r=\sqrt{17}$の方が半径4の場合に比べて二乗誤差が小さくなっています．
この理由は数式を見れば明らかです．直線フィッティングで用いる最小二乗では点と直線の差の二乗誤差を計算していますが，こちらでは各点の円の中心からの距離の二乗と半径の二乗の差の二乗誤差を計算しています．なんだか変な言い回しになってしまいましたが，先に各要素を二乗してから差を取ってその二乗誤差を計算するので，必然的に中心からの距離が大きい要素に円が引っ張られることになってしまいます．
ばらつきの少ない点群では十分な精度でフィッティングをしてくれるので，実用上はあまり気にすることはないと思いますが，通常の(感覚に沿う)最小二乗法に比べて円の外側の外れ値から受ける影響が大きいことを認識しておきましょう．
円周からの距離の誤差をもって最小二乗を計算したい場合は以下の式を最小化する必要があります．
\sum^{n}_{i=1}(\sqrt{(x_{i}-a)^{2}+(y_{i}-b)^{2}}-r)^{2} ⇒ min
別解(+ロバスト推定)
別解を見つけたので一緒に紹介したいと思います．後述しますが，こちらでは先ほどと同様の解を求められるうえ重み付きでの回帰も計算できます．(が，本記事では重みの計算までは行わず先ほどと同様の結果が得られることを確認するに留めます．)
近似計算（ロバスト推定）のアルゴリズム - 株式会社アイディール
⇒ P.6 4.2 円のロバスト推定
計算手順はまたリンク先を参考にしていただくとして，最後の計算結果
\boldsymbol{x}=(\boldsymbol{A}^{T}\boldsymbol{WA})^{-1}\boldsymbol{A}^{T}\boldsymbol{Wv}
の$\boldsymbol{W}=\boldsymbol{I}$(単位行列)として解くことで円のパラメータ$x_{0},y_{0},r$が求まります．
これで解決！と言いたいところですが，$(\boldsymbol{A}^{T}\boldsymbol{A})^{-1}$は点群のデータ数のRankを持ちます．データ数が膨大になりえることを考えれば，この逆行列をそのまま$O(n^{3})$の計算で解くのは効率的ではありません．そこで行列を分解し数式を簡略化することで計算量を抑えることを考えます．
行列の分解手法はいくつかありますが，ここでは一般的な手法として知られるQR分解を用います．
チョットワカル【QR分解】
$\boldsymbol{A}=\boldsymbol{QR}$と分解することで$\boldsymbol{A}^{T}\boldsymbol{A}=\boldsymbol{R}^{T}\boldsymbol{Q}^{T}\boldsymbol{QR}$ ここでQR分解の性質から$\boldsymbol{Q}$はユニタリ行列であるため$\boldsymbol{A}^{T}\boldsymbol{A}=\boldsymbol{R}^{T}\boldsymbol{R}$となります．
先ほどの式に代入し，式変形することで
\boldsymbol{Rx}=\boldsymbol{Q}^{T}\boldsymbol{v}
でパラメーラを計算することが出来ます．
最後に動作確認をしてみましょう．せっかくなので最初に登場した(3,0),(5,0),(0,3),(0,5),(-3,0),(-5,0),(0,-3),(0,-5)の点群で計算してみます．
A = np.array([[ 3, 0, 1],
[ 5, 0, 1],
[ 0, 3, 1],
[ 0, 5, 1],
[-3, 0, 1],
[-5, 0, 1],
[ 0, -3, 1],
[ 0, -5, 1],])
v = np.sum(np.array(np.square(A.T[:2].T)), axis=1)
#QR分解
Q, R = np.linalg.qr(A)
# x = ( 2*x_0, 2*y_0, r^2-x0^2-y0^2 )
x = np.linalg.inv(R)@Q.T@v
#解
x_0 = x[0]/2 #円の中心座標
y_0 = x[1]/2 #円の中心座標
squared_r = x[2]-x_0**2-y_0**2
print("推定された円の半径の二乗:", squared_r)
###結果###
推定された円の半径の二乗: 17.000000000000004
先ほどと同様に半径$r=\sqrt{17}$となることが確認されました．
ところで，途中$\boldsymbol{W}=\boldsymbol{I}$と仮定しましたが，この$\boldsymbol{W}$を各要素が重みとなる対角行列とすることで重み付き最小二乗近似が可能です．興味がありましたら，ぜひ計算してみてください．
最後に
あまり他者の記事をコピペしたくなかったので省略が多くなってしまいました．
そのせいで行ったり来たりして読みにくい記事になってしまっていたら申し訳ありません．
重みの計算についていつか追記するかも・・・？