Linear regression for multi-class classification 
*

*Linear regression can be used for binary classification where it competes with logistic regression. While the fitted values from linear regression are not restricted to lie between 0 and 1, unlike those from logistic regression that are interpreted as class probabilities, linear regression can still successfully assign class labels based on some threshold on fitted values (e.g. a threshold of 0.5).

*Logistic regression can be used for multi-class classification by applying it repeatedly as one-against-the rest classification. 


Can linear regression be used for multi-class classification in the same way (one against all, repeat for each class)?
Or is there some inherent feature that would to make it fail?
 A: Overview
I don't think that solving classification problems using linear regression is usually the best approach (see notes below), but it can be done. For multiclass problems, multinomial logistic regression would typically be used rather than a combination of multiple regular logistic regression models. By analogy, one could instead use least squares linear regression with multiple outputs.
Approach
Suppose we have training data $\big\{ (x_i,y_i) \big\}_{i=1}^n$ where each $x_i \in \mathbb{R}^d$ is an input point with class label $y_i$. Say there are $k$ classes. We can represent each label as a binary vector $y_i \in \{0,1\}^{k}$, whose $j$th entry is $1$ if point $i$ is a member of class $j$, otherwise $0$. The regression problem is to predict the vector-valued class labels as a linear function of the inputs, such that the squared error is minimized:
$$\min_W \ \sum_{i=1}^n \|y_i - W x_i\|^2$$
where $W \in \mathbb{R}^{k \times d}$ is a weight matrix and $\|\cdot\|^2$ is the squared $\ell_2$ norm. The inputs should contain a constant feature (i.e. one element of $x_i$ should always be $1$), so we don't have to worry about extra bias/intercept terms.
To predict the class for a new input $x$, compute the vector $a = W x$, where $a_i$ is the projection of the input onto the the $i$th row of $W$ (the weights for the $i$th class). Then, some rule can be applied to map the projections to a single class. For example, we could choose the class with the maximal projection: $\arg \max_i a_i$. This is loosely analogous to selecting the most probable class in multinomial logistic regression.
Example
Here's a plot of the decision boundaries learned from a set of 2d points, using the above method. Colors represent true class labels.

Notes
This method sacrifices the principled, probabilistic approach used in multinomial logistic regression. The squared error is also an odd choice for classification problems, where we're predicting binary values (or binary vectors, as above). The issue is that the squared error penalizes large outputs, even when these ought to be considered correct. For example, suppose the true class label is $[1,0,0]$. Outputting $[2,0,0]$ (which should correspond to high confidence in the correct class) is just as costly as outputting $[0,0,1]$ (which corresponds to high confidence in the wrong class). Even if one is willing to abandon probabilistic models, there are other loss functions designed specifically for classification, like the hinge loss used in support vector machines. The main benefit of the squared error is computational efficiency. But, this doesn't seem particularly necessary in most cases, given that we can routinely solve much more complicated problems involving massive datasets. Nevertheless, one does sometimes see the squared error used in the literature for classification problems (apparently with success). Least squares support vector machines are the most prominent example that comes to mind.
Code
Matlab code to generate the example plot above. Matrices are transposed relative to the text above, since points and labels are stored as rows.
%% generate toy dataset

% how many points and classes
n = 300;
k = 3;

% randomly choose class labels (integers from 1 to k)
c = randi(k, n, 1);

% convert labels to binary indicator vectors
% Y(i,j) = 1 if point i in class j, else 0
Y = full(sparse((1:n)', c, 1));

% mean of input points in each class
mu = [
    0, 0;
    4, 0;
    0, 4
];

% sample 2d input points from gaussian distributions
% w/ class-specific means
X = randn(n, 2) + mu(c, :);

% add a column of ones
X = [X, ones(n,1)];


%% fit weights using least squares
W = X \ Y;


%% out-of-sample prediction

% generate new test points on a grid covering the training points
[xtest2, xtest1] = ndgrid( ...
    linspace(min(X(:,2)), max(X(:,2)), 501), ...
    linspace(min(X(:,1)), max(X(:,1)), 501) ...
);
X_test = [xtest1(:), xtest2(:)];

% add a column of ones
X_test = [X_test, ones(size(X_test,1), 1)];

% project test points onto weights
A_test = X_test * W;

% predict class for each test point
% choose class w/ maximal projection
[~, c_test] = max(A_test, [], 2);


%% plot

% plot decision boundary
% using contour plot of predicted class labels at grid points
figure;
contour(xtest1, xtest2, reshape(c_test, size(xtest1)), 'color', 'k');

% plot training data colored by true class label
hold on;
scatter(X(:,1), X(:,2), [], c, 'filled');

A: I think using linear regression for multi-classes classify is totally doable.
First, let's recall that the backbone of logistic regression itself is linear regression. Following is the logistic function, but it is really just a linear function wrapped sigmoid function. When we train logistic regression, we are really just tweaking the parameter of the linear function inside the sogmoid function.

Here I would like to demonstrate with sklearn in python:
import numpy as np
from sklearn.datasets import load_iris
from sklearn.linear_model import LogisticRegression, LinearRegression
from sklearn.model_selection import train_test_split
X, y = load_iris(return_X_y=True)
X_train,X_test, y_train, y_test = train_test_split(X,y,test_size=0.2)


## logistic regression
clf_logistic = LogisticRegression(random_state=0,max_iter=1000).fit(X_train, y_train)
#print(f"score of logistic method: {clf_logistic.score(X_test, y_test)}")
### the following lines do the exact same thing
y0 = np.argmax(X_test.dot(clf_logistic.coef_.T)+clf_logistic.intercept_, axis=-1)
y1 = clf_logistic.predict(X_test)
assert np.array_equal(y0,y1) == True
acc_logistic = np.sum(y0 == y_test)/y_test.shape[0]
print(f"accuracy of logistic method: {acc_logistic}")

## linear regression
clf_linear = LinearRegression().fit(X_train, y_train)
#print(f"score of logistic method: {clf_linear.score(X_test, y_test)}")
y2 = np.round(X_test.dot(clf_linear.coef_.T)+clf_linear.intercept_)
y3 = np.round(clf_linear.predict(X_test))
assert np.array_equal(y2,y3) == True
acc_linear = np.sum(y2 == y_test)/y_test.shape[0]
print(f"accuracy of linear method: {acc_linear}")


output:
accuracy of logistic method: 0.9333
accuracy of linear method: 1.0

As to which one is better, I don't know the exact answer, it may depend on the datasets.
comments:
What the python code above does is using logistic and linear regression to classify iris into 3 categories (0,1,2). Each iris data has 4 features and 1 label. Logistic regression used one vs rest to train the data and predict the result, you can look into the shape of  clf_logistic.coef_, and clf_logistic.intercept_. However, linear regression doesn't use one vs rest (you may aslo look at the shape of clf_linear.coef_, clf_linear.intercept_), it directly predicts a numerical value, simply round float point prediction result to int, you will get the predicted label. It works well, as you can see from the printed accuracy. My understanding is that even though linear regression doesn't use one vs rest, but it works in a similar way to a single neuron like  wx + b
