Why does linear non-logistic regression work as a linear classifier? What classification error does it minimize? Suppose the data has two attributes and a label -1 or 1. So, we have a three-column matrix $X$ (two attributes and a column of ones for convenience of working with matrix notation) and a column vector $y$. We can use the following formula to compute the coefficients of linear regression that minimizes the mean squared error: $b=(X^TX)^{-1}X^Ty$. 
I generated some random training data for a classification problem and computed $b$ using the above formula. Here is the linear separator based on $b$:

Why does linear regression in three dimensions result in a good linear separator for two dimensions? Does this separator minimize some classification error (and how can this be shown)?
P.S.1.: I tried to search for an answer, but came up with mostly the opposite of what I was looking for. Namely, people say that simple (i.e. non-logistic) regression does not work well for classification. In my example it does. This inconsistency is also something I'd like to understand.
P.S.2.:
Here is the code (exported from Jupyter Notebook):
# In[1]:    
import numpy as np
import matplotlib.pyplot as plt
import math

# In[2]:
def sign(x):
    res = np.sign(x)
    return np.where(res == 0, 1, res)

# The returned X is a two-column matrix
# The returned y is a one-column matrix of labels
def data():
    X = np.random.rand(200, 2)
    X[:,0:1] *= 500
    X[:,1:2] *= 200
    y = [sign(el) for el in (X @ [0.25, 1] - 150) ]
    return X, y

# In[3]:    
def extendX(X):
    Xext = np.ones((X.shape[0], 3)) # Exercise: do it generically
    Xext[:, 0:2] = X
    return Xext

# Xext -- X extended by a column of ones
# Computes the column of coefficients b
def regression(Xext, y):
    XT = np.transpose(Xext)
    return np.linalg.inv(XT @ Xext) @ XT @ y

def nMisclassified(Xext, y, b):
    return sum(sign(Xext @ b) != sign(y))    

# In[4]:
X, y = data()
colors = [('blue' if el == 1 else 'green') for el in y]
plt.scatter(X[:,0],X[:,1],c=colors)

Xext = extendX(X)
b = regression(Xext, y)
xs = np.arange(500)
ys = (xs * b[0] + b[2]) / -b[1]
plt.plot(xs, ys, color = 'red', linewidth = 4)
plt.show()

print("Misclassified: %d" % (nMisclassified(Xext, y, b)))

 A: The linear discriminant function of Fisher takes the form $a'x$ where $x$ are the variables used for discrimination. It has the property that it maximally separates the centroids, that is, it maximizes
$$  \frac{|a'\mu_1 - a'\mu_2|^2}{\textrm{Var}(a'x)} $$
As a geometric interpretation, if you project cases on the line generated by $a$ the projections would have the least overlap.
In order to perform discrimination, you have to set a theshold, i.e. a $k$ such that you assign to group 1 if $a'x < k$ and to group 2 otherwise. This is the same as setting a separating hyperplane normal to $a$, but the exact position of that separating hyperplane is determined by the value of $k$. 
What the method (and the regression analogue) give you is the orientation of the separating hyperplane, not its exact position, which you have to set: usually, taking into consideration the gravity of the two missclassifying errors (a case of group 1 in group 2 and viceversa).
Further, linear discriminant analysis assumes equal covariance matrix of the observations in the two groups: if that fails, you would be looking at a quadratic frontier rather than a linear one. 
I think my answer was right (although too brief, but see the full details in Lachenbruch's book) if the question was: "Does an ordinary linear regression give similar results than a classical linear discriminator in the case of two groups with common dispersion matrix?"
The answer was yes. This does not imply that linear discrimination (or linear regression) are any good at discriminating, particularly if the underlying assumptions are not met!
A: The premise of the question that linear regression does a good job as a classifier is wrong. Namely, when the data is less uniformly distributed, regression does not work that well. In the following counterexample, we consider two different data sets. The blue samples are kept the same across the data sets, while the green samples are modified (the number of samples is kept the same). Regression performs quite poorly as a classifier for the second data set.

Here is the code:
# coding: utf-8

# In[1]:


import numpy as np
import matplotlib.pyplot as plt
import math


# In[2]:


# The given code
nPoints = 500
scaleX = 500
scaleY = 200

# The separator is separator_k * x + y - separator_cross_y = 0
separator_k = scaleY/(scaleX * 1.6)
separator_cross_y = 0.75 * scaleY

# Exercise! Add where to slides!
def sign(x):
    res = np.sign(x)
    return np.where(res == 0, 1, res)

# The returned X is a two-column matrix
# The returned y is a one-column matrix of labels
# Pass randomSeed, e.g. 0, to always get the same data.
def dataUniform(randomSeed = None):
    np.random.seed(randomSeed)
    X = np.random.rand(nPoints, 2)
    X[:,0:1] *= scaleX
    X[:,1:2] *= scaleY

    y = sign(X @ [separator_k, 1] - separator_cross_y)
    return X, y

# The returned X is a two-column matrix
# The returned y is a one-column matrix of labels
# Pass randomSeed, e.g. 0, to always get the same data.
def dataNonUniform(randomSeed = None):
    X1, y = dataUniform(randomSeed)

    X2 = np.zeros((nPoints, 2));
    X2[:,0] = np.linspace(0, scaleX, nPoints)
    X2[:,1] = separator_cross_y - separator_k * X2[:,0] - np.random.rand(nPoints) * scaleY / 10

    X = np.ones((nPoints,2))
    X = np.where(y[:,np.newaxis] == 1, X1, X2)

    return X, y


# In[3]:


# The following two functions are the subject of the assignment

def extendX(X):
    Xext = np.ones((X.shape[0], 3)) # Exercise: do it generically
    Xext[:, 0:2] = X
    return Xext

# Xext -- X extended by a column of ones
# Computes the column of coefficients b
def regression(Xext, y):
    XT = np.transpose(Xext)
    return np.linalg.inv(XT @ Xext) @ XT @ y

def nMisclassified(Xext, y, b):
    return sum(sign(Xext @ b) != sign(y))    


# In[4]:


for gen in [dataUniform, dataNonUniform]:
    X, y = gen(0)
    colors = [('blue' if el == 1 else 'green') for el in y]
    plt.scatter(X[:,0],X[:,1],c=colors)

    Xext = extendX(X)
    b = regression(Xext, y)
    xs = np.arange(500)
    ys = (xs * b[0] + b[2]) / -b[1]
    plt.plot(xs, ys, color = 'red', linewidth = 4)
    plt.show()

    print("Misclassified: %d" % (nMisclassified(Xext, y, b)))

