Email Topic Extraction and Clustering with Non-negative Matrix Factorization

So I was checking some email and found myself inundated with not-quite-spam. I bet that you have a lot of email that is neither totally unwanted scam emails, nor personal correspondence. You probably want to pay different levels of attention to the these mails based on their content. Wouldn’t it be great if you could automate all of your email handling?

Let’s take a look at a corpus of emails and extract topic structure with Non-negative Matrix Factorization. You can use run this code on your computer to use the email currently stored on your desktop and find patterns in your email. This post uses the mails currently stored in my Maildir. It turns out that most of my email is sent by machines!

The first step is to include dependencies and define some helper functions

%matplotlib inline
import pip
pip.main(['install', 'scipy', 'sklearn', 'matplotlib', 'pandas'])

from __future__ import print_function
from time import time

import mailbox
import pandas as pd
from sklearn.feature_extraction.text import TfidfVectorizer, CountVectorizer
from sklearn.decomposition import NMF, LatentDirichletAllocation
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')
plt.rc('font', **{'family' : 'monospace',
          'weight' : 'bold',
          'size'   : 30})

n_samples = 40000
n_features = 2000
n_topics = 10
n_top_words = 30

Accessing the Email

The Python standard library mailbox module provides easy access to emails in the Maildir format. You can use the great tool isync/mbsync to synchronize your email to your desktop. I have used it on Linux and OSX.

md = mailbox.Maildir("/home/james/Maildir/jjpfnet/INBOX")

# The code in this cell (and only this cell) is modified from @Tyndyll's 
# gist https://gist.github.com/tyndyll/6f6145f8b1e82d8b0ad8
# which is available on gist.github.com
def extract(msg, msg_obj, attachments):
    if msg.is_multipart():
        for part in msg.get_payload():
            if part.is_multipart():
                extract(part, msg_obj, attachments)
                continue
            if part.get('Content-Disposition') is None:
                msg_obj["body"] += part.get_payload()
            else:
                if part.get('Content-Disposition').startswith('attachment'):
                    attachments[part.get_filename()] = {
                            'data': part.get_payload(decode=True),
                            'mime': part.get_content_type()
                        }
    else:
        msg_obj["body"] += msg.get_payload()
    return

"""
    messages(mdir) is a generator for message objects from a maildir object.
"""
def messages(mdir, directory):
    for msg in mdir:
        msg_obj = {"path": directory, "body": ""}
        attachments = {}
        for key, value in msg.items():
            msg_obj[key.lower()] = value
        lists = ["from", "cc", "bcc", "to"]
        extract(msg, msg_obj, attachments)
        yield msg_obj, attachments
        

All you have to do is open the mailbox and pull out all the subject and body text. In order for an Machine Learning tool to work you have to filter out some noise. Email has a lot of junk that shows up in some machine generated emails. The most important thing is to remove as much of the random strings that aren’t words from the corpus. Minified HTML,JS, and CSS also make processing more difficult.

This filtering step would be much more involved in a production application of data analysis. You would have to do the heavy lifting of extracting as much useful text as possible while discarding the noise.

i = 0
bodies = []
allsubs = []
body = ''
mesgs = messages(md, "INBOX")
for mesg, attach in mesgs:
    try:
        # we need to filter out long words because they are either
        # minified javascript/css or API keys
        body = ' '.join(filter(lambda x: len(x) <25, mesg['body'].split()))
    except Exception as ex:
        print('No body found for message {}'.format(mesg))
        print('\tException: {}'.format(ex))
    bodies.append(body)
    allsubs.append(mesg['subject'])
    i+=1
    if i > n_samples:
        break

# then we throw out documents with less than 2 words
# this catches a lot of emails that were completely composed of
# minified css and js
data_samples = list(filter(lambda x: (len(x.split()) > 2), bodies))
subs = []
for i, bdy in enumerate(bodies):
    if len(bdy.split()) > 2:
        subs.append(allsubs[i])
print('The number of documents is: {}'.format(len(data_samples)))

The number of documents is: 4263

4,000 emails is not that much, but we can still apply these techniques. A good comparison would be a deep neural net, which might take 1-2 orders of magnitude more data to train.

Take a look at some of my email for fun. You should be able to eyeball the difference between mail that came from an automated system (potential spam) from email that came from humans.

This is the account that receives my github notifications so a lot of the email is about code.

for i, x in enumerate(data_samples[0:10]):
    print('Subject: {}'.format(subs[i]))
    print('Body:')
    print(x[:min(240, len(x))])
    print('------------------------------Message Truncated-----------------\n')

    Subject: Re: Pollux has an inconsistent root file system.
    Body:
    I am also forwarding this email to Anita who has root privileges on the cas= tor and pollux. Best, Oded From: Jason Riedy Sent: Monday, April 10, 2017 2:52 PM To: Ediger, David M; Robert McColl; James Fairbanks; Green, Oded Subject: Pollux 
    ------------------------------Message Truncated-----------------
    
    Subject: Purchase order -241PO1800282
    Body:
    Dear Sir, Please find attached purchase order for your reference. Kindly qoute your best prices with proforma invoice, payment terms and delivery schedule. Thanks & Best Regards Mr. Jafar Sadiq Sales Excutive Mob:  S
    ------------------------------Message Truncated-----------------
    
    Subject: Re: How is the atomic arrays going 
    Body:
    Hi James, I wrote a package implementing the functionality called UnsafeAtomics.jl<ht= I just ported all the relev= ant atomic tests from Base.Threads tests to the unsafe package and they all= pass. I think I will post on the list as a foll
    ------------------------------Message Truncated-----------------
    
    Subject: Re: [JuliaGraphs/LightGraphs.jl] should is_cyclic return true on a
     2-node PathGraph? (#592)
    Body:
    Closed #592. -- You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub: <a class="issue-link js-issue-link" data-id="221955592" data-error-text="Failed to load issue title" title is private">#59
    ------------------------------Message Truncated-----------------
    
    Subject: RE : Potential Position for Senior DevOps Engineer/System Administrator Loc: Princeton, NJ Dur: 6+ M Rate: OPEN
    Body:
    <p>Hi, </p> <p> I received your contact from a referral as a potential for contract / full time position. Based on what I learnt, I feel you might be suitable for the following position as given in the job description attached at the bottom
    ------------------------------Message Truncated-----------------
    
    Subject: Re: [jpoovey3/HDAPWorkflows] Move execution to self.run method; Add
     docstrings; (#1)
    Body:
    @jfairbanks6, did you want to go ahead and merge this one? I've been working on another branch `plp_model` for finishing up the modeling workflow. -- You are receiving this because you were mentioned. Reply to this email directly or view it
    ------------------------------Message Truncated-----------------
    
    Subject: Re: [JuliaGraphs/LightGraphs.jl] WIP: LightGraphs Abstraction (#541)
    Body:
    @sbromberger pushed 1 commit. d8d12ce 1.0 -> 0.8 -- You are receiving this because you are subscribed to this thread. View it on GitHub: <p><a pushed 1 commit.</p> <ul> <li><a 1.0 -&gt; 0.8</li> </ul> <p />You are receiving this because you
    ------------------------------Message Truncated-----------------
    
    Subject: Re: [stingergraph/stinger] compile without warnings (#218)
    Body:
    @zeuner pushed 1 commit. 3eaf610 linking related fixes -- You are receiving this because you are subscribed to this thread. View it on GitHub: <p><a pushed 1 commit.</p> <ul> <li><a linking related fixes</li> </ul> <p />You are receiving th
    ------------------------------Message Truncated-----------------
    
    Subject: Re: [JuliaGraphs/LightGraphs.jl] WIP: Eager implementation of Prim's
     algorithm (#530)
    Body:
    Oh! I hope you're feeling better. -- You are receiving this because you were mentioned. Reply to this email directly or view it on GitHub: I hope you're feeling better.</p> <p />You are receiving this because you were mentioned.<br />Reply 
    ------------------------------Message Truncated-----------------
    
    Subject: Re: [JuliaGraphs/LightGraphs.jl] merge vertices operator (#676)
    Body:
    # Report > Merging into will **decrease** coverage by `0.39%`. > The diff coverage is `0%`. ```diff @@ Coverage Diff @@ ## master #676 +/- ## - Coverage 100% 99.6% -0.4% Files 59 59 Lines 3026 3038 +12 Hits 3026 3026 - Misses 0 12 +12 ``` -
    ------------------------------Message Truncated-----------------

Running the Model

After all that work to get the data ready, it is time to fit a model to the data! Good thing our friends at Sklearn have made NMF and TF-IDF feature extraction a simple library call. TF-IDF is the matrix that is derived from a corpus of text where each entry of a matrix $A$ is defined as the $A[i,j]$ is the frequency of the $j$th term in the $i$th document divided by the number of documents that contains term $j$. This matrix captures the idea that common words are less semantically important than rare words.

This code adapted from the Sklearn documentation.

# Use tf-idf features for NMF.
print("Extracting tf-idf features for NMF...")
tfidf_vectorizer = TfidfVectorizer(max_df=0.95, min_df=2,
                                   max_features=n_features,
                                   stop_words='english')
t0 = time()
tfidf = tfidf_vectorizer.fit_transform(data_samples)
print("done in %0.3fs." % (time() - t0))

# Fit the NMF model
print("Fitting the NMF model with tf-idf features, "
      "n_samples=%d and n_features=%d..."
      % (n_samples, n_features))
t0 = time()
nmf = NMF(n_components=n_topics, random_state=1,
          alpha=.1, l1_ratio=.5).fit(tfidf)
print('NMF took {} iterations of CD'.format(nmf.n_iter_))
print("done in %0.3fs." % (time() - t0))

tfidf_feature_names = tfidf_vectorizer.get_feature_names()

Extracting tf-idf features for NMF...
done in 2.683s.
Fitting the NMF model with tf-idf features, n_samples=40000 and n_features=2000...
NMF took 134 iterations of CD
done in 0.931s.

Now we want to visualize the topics. The first step is to view them as lists of words. The machine has only learned sets of words it doesn’t have any idea of conceptual similarity of the words. They are just words that make good low rank factors to the TF-IDF matrix. Fortunately your brian will identify “topics” when you read the lists of words. Or, at least it would if this were your email!

def top_words(topic, n_top_words):
    return topic.argsort()[:-n_top_words - 1:-1]
def topic_table(model, feature_names, n_top_words):
    topics = {}
    for topic_idx, topic in enumerate(model.components_):
        t = ("topic_%d:" % topic_idx)
        topics[t] = [feature_names[i] for i in top_words(topic, n_top_words)]
    return pd.DataFrame(topics)

topic_table(nmf, tfidf_feature_names, n_top_words).head(10)

Analyzing the Results

The results of the matrix factorization can give us the list of words in each topic, but we want to know more about the structure of the data. Is there structure between the topics? How do the documents break down in terms of topics? Are there clusters of documents? Can we detect spam with this model?

gittopics = [0,2,5,7,8]
csstopics = [1,3,4,6]

docweights = nmf.transform(tfidf_vectorizer.transform(data_samples))

plt.figure(figsize=(40,10))
plt.imshow(docweights[0:200, :].T, aspect=8, interpolation='none')
plt.title('Topic Weighting from NMF')
plt.xlabel('document')
plt.ylabel('Topic')
plt.colorbar()

We can also see the association between the topics by computing the inner product of the normal equations for the topic-document matrix.

$$ A \approx HW $$

$A$ is documents by words $H$ is documents by topics and $W$ is topics by words.

The product $H’H$ is the association of the topics with each other based on occuring in the same documents.

The product $WW’$ is the association of the topics with each other based on their overlapping distributions of words.

A = tfidf_vectorizer.transform(data_samples)
W = nmf.components_
H = nmf.transform(A)

print('A is {} x {}'.format(*A.shape))
print('W is {} x {}'.format(*W.shape))
print('H is {} x {}'.format(*H.shape))

A is 4263 x 2000
W is 10 x 2000
H is 4263 x 10

Here we see that there are sets of topics that share a lot of documents. These topics are related to each other.

Topics gittopics share a lot of documents.

Topics csstopics share a lot of documents, they are topics related to code.

These associations will come up in more detail later.

print('Git Topics: ', gittopics)
print('CSS Topics: ', csstopics)

Git Topics:  [0, 2, 5, 7, 8]
CSS Topics:  [1, 3, 4, 6]

plt.figure(figsize=(20,20))
plt.imshow(H.T.dot(H), interpolation='none')
plt.colorbar()
plt.title("Association of Topics (H'H)")
plt.xlabel('Topic index')

Here you can see that the words offer some structure to the topics.

Topics 5 and 6 overlap which matches the intuition that they are both sets of CSS related words.

However topics 0,2,4, and 7 have some overlap as they are topics related to code. Topic 9 is about LLNL

plt.figure(figsize=(20,20))
plt.imshow(W.dot(W.T), interpolation='none')
plt.colorbar()
plt.title("Association of Topics (WW')")
plt.xlabel('Topic index')

Spot Checking the Data: Annecdatal evidence

One of the reasons that Machine Learning is so popular with Images, Text, and Audio is that humans are equiped for making judgements about accuracy in these media. By comparison, mechanical failures, geospatial data, medical records are all complicated domains where humans don’t have a natural ability to judge correctness. The other reason, of course, is the ubiquity of Images, Text, and Audio on the web.

So let’s use our human abilities to recognize accuracy on this clustering problem and spot check the clusters.

Topic models yield a clustering technique by associating each document with its dominant topic and then looking at representative examples from each cluster.

df = pd.DataFrame({'subject':subs,
                   'topic':docweights.argmax(axis=1),
                   'body':data_samples},
                  columns=['topic', 'subject', 'body'])

"""
    showdocs(df, topics, nshow=5) is a function that gathers a number of 
    documents from a set of topics as a dataframe.
"""
def showdocs(df, topics, nshow=5):
    idx = df.topic == topics[0]
    for i in range(1, len(topics)):
        idx = idx | (df.topic == topics[i])
    return df[idx].groupby('topic').head(nshow).sort_values('topic')

These topics are mostly from Newsletters and are identified by their CSS words like font, 0pt, width, padding, margin, quot, align, serif.

showdocs(df, csstopics)

These topics are mostly from Github and relate to code, with words like pull, commit, julia, graph, mergining, application.

showdocs(df, gittopics, nshow=3)

Analysis of Residuals

We can look at the quality of the model on various documents. For some documents the text is well approximated by the latent factors, but for some other documents the text is poorly approximated by the topics.

The residual for a document is defined as $$ r[i] = ||A[i,:] - H[i,:]W||_2 $$

When $r[i]$ is small the model is capturing the text of the document, when the residual is large, the document is using some combination of words that does not lie on the subspace spaned by the word vectors (columns of W).

We see that the residuals are small for documents that fit the mental category that we assign to the topic and large for documents that are individualized and personal.

import numpy as np
r = np.zeros(A.shape[0])

for row in range(A.shape[0]):
    r[row] = np.linalg.norm(A[row, :]- H[row,:].dot(W))

print(sum(np.sqrt(r)))
nmf.reconstruction_err_

3615.85321153
50.592531594804974

plt.figure(figsize=(20,20))
plt.hist(r);
plt.title('Histogram of residuals')
plt.xlabel('$r_i$')
plt.ylabel('count($i$ | $r[i] = x$)')

df['residual'] = r
df.sort_values(['topic','residual'],
                             ascending=True).groupby('topic').head(5)

I would put some large residual examples in this presentation, but they are mostly non-spam and have been excluded for privacy. You can run this experiment on your own emails and find that you can build a pretty good spam filter and message classifier from topic modeling with non-negative matrix factorization.

Conclusion

We can see that email inboxes (at least mine) contain several distinct topics of messages, these messages can be detected Non-negative Matrix Factorization. There is reason to believe that the residual of a topic model can be used to build a spam filter. While verification of the utility of NMF topic modeling for spam detection requires labeled data, spot checks suggest this will work.

One can learn about the structure of the topics and documents by analyzing the factors $A \approx H,W$ to gain insight into the salient features of the data detected by the algorithm.