Compute and graph the LDA decision boundary I saw an LDA (linear discriminant analysis) plot with decision boundaries from The Elements of Statistical Learning: 
I understand that data are projected onto a lower-dimensional subspace. However, I would like to know how we get the decision boundaries in the original dimension such that I can project the decision boundaries onto a lower-dimensional subspace (likes the black lines in the image above).
Is there a formula that I can use to compute the decision boundaries in the original (higher) dimension? If yes, then what inputs does this formula need?
 A: This particular figure in Hastie et al. was produced without computing equations of class boundaries. Instead, algorithm outlined by @ttnphns in the comments was used, see footnote 2 in section 4.3, page 110:

For this figure and many similar figures in the book we compute the decision boundaries by an exhaustive contouring method. We compute the decision rule on a fine lattice of points, and then use contouring algorithms to compute the boundaries.

However, I will proceed with describing how to obtain equations of LDA class boundaries.
Let us start with a simple 2D example. Here is the data from the Iris dataset; I discard petal measurements and only consider sepal length and sepal width. Three classes are marked with red, green and blue colours:

Let us denote class means (centroids) as $\boldsymbol\mu_1, \boldsymbol\mu_2, \boldsymbol\mu_3$. LDA assumes that all classes have the same within-class covariance; given the data, this shared covariance matrix is estimated (up to the scaling) as $\mathbf{W} = \sum_i (\mathbf{x}_i-\boldsymbol \mu_k)(\mathbf{x}_i-\boldsymbol \mu_k)^\top$, where the sum is over all data points and centroid of the respective class is subtracted from each point.
For each pair of classes (e.g. class $1$ and $2$) there is a class boundary between them. It is obvious that the boundary has to pass through the middle-point between the two class centroids $(\boldsymbol \mu_{1} + \boldsymbol \mu_{2})/2$. One of the central LDA results is that this boundary is a straight line orthogonal to $\mathbf{W}^{-1} \boldsymbol (\boldsymbol \mu_{1} - \boldsymbol \mu_{2})$. There are several ways to obtain this result, and even though it was not part of the question, I will briefly hint at three of them in the Appendix below.
Note that what is written above is already a precise specification of the boundary. If one wants to have a line equation in the standard form $y=ax+b$, then coefficients $a$ and $b$ can be computed and will be given by some messy formulas. I can hardly imagine a situation when this would be needed.
Let us now apply this formula to the Iris example. For each pair of classes I find a middle point and plot a line perpendicular to $\mathbf{W}^{-1} 
\boldsymbol (\boldsymbol \mu_{i} - \boldsymbol \mu_{j})$:

Three lines intersect in one point, as should have been expected. Decision boundaries are given by rays starting from the intersection point:

Note that if the number of classes is $K\gg 2$, then there will be $K(K-1)/2$ pairs of classes and so a lot of lines, all intersecting in a tangled mess. To draw a nice picture like the one from the Hastie et al., one needs to keep only the necessary segments, and it is a separate algorithmic problem in itself (not related to LDA in any way, because one does not need it to do the classification; to classify a point, either check the Mahalanobis distance to each class and choose the one with the lowest distance, or use a series or pairwise LDAs).
In $D>2$ dimensions the formula stays exactly the same: boundary is orthogonal to $\mathbf{W}^{-1} \boldsymbol (\boldsymbol \mu_{1} - \boldsymbol \mu_{2})$ and passes through $(\boldsymbol \mu_{1} + \boldsymbol \mu_{2})/2$. However, in higher dimensions this is not a line anymore, but a hyperplane of $D-1$ dimensions. For illustration purposes, one can simply project the dataset to the first two discriminant axes, and thus reduce the problem to the 2D case (that I believe is what Hastie et al. did to produce that figure).
Appendix
How to see that the boundary is a straight line orthogonal to $\mathbf{W}^{-1} (\boldsymbol \mu_{1} - \boldsymbol \mu_{2})$? Here are several possible ways to obtain this result:


*

*The fancy way: $\mathbf{W}^{-1}$ induces Mahalanobis metric on the plane;  the boundary has to be orthogonal to $\boldsymbol \mu_{1} - \boldsymbol \mu_{2}$ in this metric, QED.

*The standard Gaussian way: if both classes are described by Gaussian distributions, then the log-likelihood that a point $\mathbf x$ belongs to class $k$ is proportional to $(\mathbf x - \boldsymbol \mu_k)^\top \mathbf W^{-1}(\mathbf x - \boldsymbol \mu_k)$. On the boundary the likelihoods of belonging to classes $1$ and $2$ are equal; write it down, simplify, and you will immediately get to $\mathbf x^\top \mathbf W^{-1} (\boldsymbol \mu_{1} - \boldsymbol \mu_{2}) = \mathrm{const}$, QED.

*The laboursome but intuitive way. Imagine that $\mathbf{W}$ is an identity matrix, i.e. all classes are spherical. Then the solution is obvious: boundary is simply orthogonal to $\boldsymbol \mu_1 - \boldsymbol \mu_2$. If classes are not spherical, then one can make them such by sphering. If the eigen-decomposition of $\mathbf{W}$ is $\mathbf{W} = \mathbf U \mathbf D \mathbf U^\top$, then matrix $\mathbf S = \mathbf D^{-1/2} \mathbf U^\top$ will do the trick (see e.g. here). So after applying $\mathbf S$, the boundary is orthogonal to $\mathbf S (\boldsymbol \mu_{1} - \boldsymbol \mu_{2})$. If we take this boundary, transform it back with $\mathbf S^{-1}$ and ask what is it now orthogonal to, the answer (left as an exercise) is: to $\mathbf S^\top \mathbf S \boldsymbol (\boldsymbol \mu_{1} - \boldsymbol \mu_{2})$. Plugging in the expression for $\mathbf S$, we get QED.
A: I want to start of by thanking @amoeba says Reinstate Monica & ttnphns for their contributions that have greatly helped me! I'm so grateful in fact, I'd want to buy them a drink or be able to return the favor somehow.
The only thing I'm going to add to their reply is my python implementation of drawing these Decision boundaries, I think it'll help others, theory and insight is great, but some understand better through code.
The code below is useful for visualization, I have used LDA for dimensionality reduction (10 000 dim to 2D) for 3 classes. The framework is sklearn. Only the code to plot the DB is written below, if you're interested in the training part of the classifier, sci-kit learn's documentation is VERY good. The code below assumes you have projected your training data to 2D and you have computed their means.
x_min, x_max = plt.xlim()
# xplot = np.linspace(x_min, x_max, 100) #to plot the whole lines and find the intersection
#x_l1 and x_l2 to only plot the relevant part of the lines [found the intersection by first plotting the data]
x_l1 = np.linspace(x_min, 0.272, 100)
x_l2 = np.linspace(0.272, x_max, 100)

cov = lda1_2features.covariance_
prec_m = np.linalg.inv(cov)    

line1 = np.dot(prec_m, (mu1[0]-mu1[1]).T)#mu1 contains the coordinates of all the
# classes [in this example 3 classes)
midpoint1 = (mu1[1]+mu1[0])/2
plt.plot(midpoint1[0], midpoint1[1],
        'p', color='magenta', markersize=10, markeredgecolor='grey')
rico1 = -line1[0]/line1[1]
cte1 = midpoint1[1]-(rico1)*midpoint1[0]
# plt.plot(xplot, (rico1*xplot)+cte1, '--b')
plt.plot(x_l1, (rico1*x_l1)+cte1, '--b')

line2 = np.dot(prec_m, (mu1[0]-mu1[2]).T)
midpoint2 = (mu1[2]+mu1[0])/2
plt.plot(midpoint2[0], midpoint2[1],
        'p', color='magenta', markersize=10, markeredgecolor='grey')
rico2 = -line2[0]/line2[1]
cte2 = midpoint2[1]-(rico2)*midpoint2[0]
# plt.plot(xplot, (rico2*xplot)+cte2, '--r')
plt.plot(x_l2, (rico2*x_l2)+cte2, '--r')

line3 = np.dot(prec_m, (mu1[1]-mu1[2]).T)
midpoint3 = (mu1[2]+mu1[1])/2
plt.plot(midpoint3[0], midpoint3[1],
        'p', color='magenta', markersize=10, markeredgecolor='grey')
rico3 = -line3[0]/line3[1]
cte3 = midpoint3[1]-(rico3)*midpoint3[0]
# plt.plot(xplot, (rico3*xplot)+cte3, '--g')
plt.plot(x_l1, (rico3*x_l1)+cte3, '--g')
plt.show()


2D subspace of classifier I with Decision boundaries, class centers(yellow star), cov_ellispses and midpoints(pentagons) plotted in one figure using matplotlib]
Sources;
scikit_lda_qda
scikit_lda
