The green line 'a' represents my original
approximated best fit line using sliders m_1 and b_1.
The red line 'k' represents the calculated best fit
line that minimizes the squares of the distances
of the four points to the line.
Notice that when the greeen line is brought as close
to the red line, the square of the distances in the 
spreadsheet is at a minimum.

The pink line 'l' represents the best fit line when vertical
distances of the points to the line are minimized.
When all four points are aligned vertically, both lines
become undefined because their calculated values
end up dividing by zero.

A = (-10.762614546822439, -3.785040299743387)
B = (4.848775761083326, 3.245621292824141)
C = (-8.326794782468347, -2.5671304175663345)
D = (3.686225419005237, 6.622553238860513)
m = 0.5 (slider from -5 to 5, step 0.1)
b = 1.5 (slider from -5 to 5, step 0.1)
a: y = m * x + b
c = OrthogonalLine(B, yAxis)
d = OrthogonalLine(A, yAxis)
e = OrthogonalLine(C, yAxis)
f = OrthogonalLine(D, yAxis)
G = Intersect(a, e)
i = Segment(C, G)
H = Intersect(a, f)
j = Segment(D, H)
E = Intersect(a, d)
g = Segment(A, E)
F = Intersect(a, c)
h = Segment(B, F)
A1: x(A)
A2: x(B)
A3: x(C)
A4: x(D)
B1: y(A)
B2: y(B)
B3: y(C)
B4: y(D)
C1: x(A) - x(E)
C2: x(B) - x(F)
C3: x(C) - x(G)
C4: x(D) - x(H)
D1: C1^(2)
D2: C2^(2)
D3: C3^(2)
D4: C4^(2)
D6 = Sum(D1:D4)
m_h: (-3 * y(A)^(2) + 2 * y(A) * y(B) - 3 * y(B)^(2) + 2 * y(A) * y(C) + 2 * y(B) * y(C) - 3 * y(C)^(2) + 2 * y(A) * y(D) + 2 * y(B) * y(D) + 2 * y(C) * y(D) - 3 * y(D)^(2)) / (-3 * x(A) * y(A) + x(B) * y(A) + x(C) * y(A) + x(D) * y(A) + x(A) * y(B) - 3 * x(B) * y(B) + x(C) * y(B) + x(D) * y(B) + x(A) * y(C) + x(B) * y(C) - 3 * x(C) * y(C) + x(D) * y(C) + x(A) * y(D) + x(B) * y(D) + x(C) * y(D) - 3 * x(D) * y(D))
b_h: (-x(B) * y(A)^(2) - x(C) * y(A)^(2) - x(D) * y(A)^(2) + x(A) * y(A) * y(B) + x(B) * y(A) * y(B) - x(A) * y(B)^(2) - x(C) * y(B)^(2) - x(D) * y(B)^(2) + x(A) * y(A) * y(C) + x(C) * y(A) * y(C) + x(B) * y(B) * y(C) + x(C) * y(B) * y(C) - x(A) * y(C)^(2) - x(B) * y(C)^(2) - x(D) * y(C)^(2) + x(A) * y(A) * y(D) + x(D) * y(A) * y(D) + x(B) * y(B) * y(D) + x(D) * y(B) * y(D) + x(C) * y(C) * y(D) + x(D) * y(C) * y(D) - x(A) * y(D)^(2) - x(B) * y(D)^(2) - x(C) * y(D)^(2)) / (3 * x(A) * y(A) - x(B) * y(A) - x(C) * y(A) - x(D) * y(A) - x(A) * y(B) + 3 * x(B) * y(B) - x(C) * y(B) - x(D) * y(B) - x(A) * y(C) - x(B) * y(C) + 3 * x(C) * y(C) - x(D) * y(C) - x(A) * y(D) - x(B) * y(D) - x(C) * y(D) + 3 * x(D) * y(D))
k: y = m_h * x + b_h
m_v: (-(-(2 * x(A) + 2 * x(B) + 2 * x(C) + 2 * x(D)) * (-2 * y(A) - 2 * y(B) - 2 * y(C) - 2 * y(D)) + 8 * (-2 * x(A) * y(A) - 2 * x(B) * y(B) - 2 * x(C) * y(C) - 2 * x(D) * y(D)))) / (-(2 * x(A) + 2 * x(B) + 2 * x(C) + 2 * x(D))^(2) + 8 * (2 * x(A)^(2) + 2 * x(B)^(2) + 2 * x(C)^(2) + 2 * x(D)^(2)))
b_v: 1 / 4 * (y(A) + y(B) + y(C) + y(D)) + ((x(A) / 4 + x(B) / 4 + x(C) / 4 + x(D) / 4) * (-(2 * x(A) + 2 * x(B) + 2 * x(C) + 2 * x(D)) * (-2 * y(A) - 2 * y(B) - 2 * y(C) - 2 * y(D)) + 8 * (-2 * x(A) * y(A) - 2 * x(B) * y(B) - 2 * x(C) * y(C) - 2 * x(D) * y(D)))) / (-(2 * x(A) + 2 * x(B) + 2 * x(C) + 2 * x(D))^(2) + 8 * (2 * x(A)^(2) + 2 * x(B)^(2) + 2 * x(C)^(2) + 2 * x(D)^(2)))
l: y = m_v * x + b_v