# 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?

• Rather than decision boundaries you will probably find more utility in considering posterior probabilities of class membership. This can be done with fewer assumptions using polytomous (multinomial) logistic regression but can also be done with LDA (posterior probabilities). Apr 2, 2014 at 12:28
• Within LDA, those classification boundaries constitute what is known a territorial map. I work with SPSS, and it plots it, although in text format. According to one SPSS designer, the boundaries are found easily by practical approach: Jun 18, 2014 at 9:37
• (cont.) every point of a fine grid is LDA-classified, and then if a point was classified as its neighbours were, that point is not shown. Thus only boundaries as "bands of ambiguety" are left in the end. Citation: they (bondaries) are never computed. The plot is drawn by classifying every character cell in it, then blanking out all those surrounded by cells classified into the same category. Jun 18, 2014 at 9:41

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:

1. 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.

2. 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.

3. 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.

• I have not been studying your answer. It seems sophisticated and may be correct. What is about the practical and easier "sprinkle points, classify, then deduce boundaries" approach that I outlined in a comment? Is your approach comparable with its results (which are obviously correct)? What do you think? Jun 18, 2014 at 9:51
• @ttnphns: The only technical part of my answer (a numbered list with 3 items) is providing some proofs and can safely be skipped. The rest, I believe, is not particularly sophisticated! Maybe I should move that "extra" part down, as an appendix? Regarding your comments: I think this is a valid approach, and I like the ASCII looks of the SPSS "territorial map". Maybe you could move your comments into a separate answer (and give an exemplary picture of the SPSS map there), I think it would be helpful for future references. The results should of course be equivalent. Jun 18, 2014 at 21:48
• @ttnphns: It turns out that Hastie et al. used exactly the method you described here to plot their figures, including the one reproduced in the OP. I found a footnote saying exactly that (and updated my answer, quoting it in the beginning). Mar 16, 2015 at 22:54
• Waouh! excellent answer (3 years later !) may I ask how you got to draw the segments in this particular problem ? Jun 15, 2018 at 22:58
• Great answer! Let me inform you that in order for the statement "The boundary has to pass through the middle-point between the two class centroids $(\mu_1 + \mu_2)/2$." to be true, we need the additional assumption $\pi_1 = \pi_2$ P.S. Mar 6, 2020 at 5:24

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