《集体智慧编程》读书笔记9

最近重读《集体智慧编程》，这本当年出版的介绍推荐系统的书，在当时看来很引领潮流，放眼现在已经成了各互联网公司必备的技术。
这次边阅读边尝试将书中的一些Python语言例子用C#来实现，利于自己理解，代码贴在文中方便各位园友学习。

由于本文可能涉及到的与原书版权问题，请第三方不要以任何形式转载，谢谢合作。

第九部分

这一部分介绍是一种非监督学习方法，作用可总结为在一个数据集中提取出重要的潜在特征。这项工作又称为特征提取。特征提取的使用场景包括在混杂的声音中分离出单独的声音，以及在一组文档中识别出单词的组合形式从而在文档中识别出主题。

提取新闻主题

作为本部分的示例，我们来详细看看如何在一组文章中去提取主题，这些文章来自不同的订阅源。
最终可以看到这些文章中哪些有共同的主题，或者一篇文章可能包含不同的主题。

获取新闻文章

我们通过书中提供的那些Rss源来获取文章，本书出版到现在至少已经7年了。这其中有一部分Rss源已经变化。但大部分经过一些调整依旧可以访问。楼主诸个尝试了下，把依旧存在的源整理了出来。

如果园友遇到不能访问的链接说明你可能需要一个梯子。

先项目中添加一个名为NewsFeatures的类，并将Rss源的链接作为静态成员放入其中。

public class NewsFeatures
{
    private static readonly List<string> Feedlist = new List<string>()
    {
        "http://feeds.reuters.com/reuters/topNews",
        "http://feeds.reuters.com/Reuters/domesticNews",
        "http://feeds.reuters.com/Reuters/worldNews",
        "http://rss.nytimes.com/services/xml/rss/nyt/HomePage.xml",
        "http://rss.nytimes.com/services/xml/rss/nyt/International.xml",
        "https://news.google.com/?output=rss",
        "http://www.foxnews.com/xmlfeed/rss/0,4313,0,00.rss",
        "http://rss.cnn.com/rss/edition.rss",
        "http://rss.cnn.com/rss/edition_world.rss",
        "http://rss.cnn.com/rss/edition_us.rss"
    };
}

.NET平台最常用的RSS解析工具就是.NET Framework3.5开始System.ServiceModel内置的SyndicationFeed，但是这个类的一个问题就是容错性不佳。处理上面的那些RSS源就经常遇到xml中时间等无法解析而报错的情况。
所以这里使用了一个名为FeedReader的第三方库来完成RSS的解析。

public class FeedParser
{
    public static Feed Parse(string url)
    {
        try
        {
            // 这种反异步的方式应该只用于测试中
            var feed = FeedReader.ReadAsync(url).GetAwaiter().GetResult();
            return feed;
        }
        catch (Exception e)
        {
            return null;
        }
    }
}

为了清洗由RSS得到的内容，需要一个函数来去除Html标签，将其加入到NewsFeatures中。

public string StripHtml(string h)
{
    var psb = new StringBuilder();
    var s = 0;
    foreach (var c in h)
    {
        if (c == '<')
            s = 1;
        else if (c == '>')
        {
            s = 0;
            psb.Append(" ");
        }
        else if (s == 0)
            psb.Append(c);
    }
    return psb.ToString();
}

然后我们剔除一些无用字符，长度过短的单词。经过这个清洗后，从网络下载的文章就基本可以用于接下来的分析。

public List<string> SeparateWords(string text)
{
    return Regex.Split(text, "\W")
        .Where(s => s.Length > 3)
        .Select(s => s.ToLower()).ToList();
}

将上面的方法组合一下，从订阅源下载新闻，并使用上面的函数对内容进行清洗。

public (Dictionary<string, int>, Dictionary<int, Dictionary<string, int>>,
    List<string>) GetArticleWords()
{
    var allWords = new Dictionary<string, int>();
    var articleWords = new Dictionary<int, Dictionary<string, int>>();
    var articleTitles = new List<string>();
    var ec = 0;
    // 遍历所有RssFeed
    foreach (var feed in Feedlist)
    {
        var f = FeedParser.Parse(feed);

        // 遍历每篇文章
        foreach (var e in f.Items)
        {
            //跳过标题相同的文章
            if (articleTitles.Contains(e.Title))
                continue;
            // 提取单词
            var txt = e.Title + StripHtml(e.Description??string.Empty);
            var words = SeparateWords(txt);
            articleTitles.Add(e.Title);

            // 在allWords和articleWords中增加针对当前单词的计数
            foreach (var word in words)
            {
                if (!allWords.ContainsKey(word))
                    allWords.Add(word, 0);
                allWords[word] += 1;
                if (!articleWords.ContainsKey(ec))
                    articleWords.Add(ec, new Dictionary<string, int>());
                if (!articleWords[ec].ContainsKey(word))
                    articleWords[ec].Add(word, 0);
                articleWords[ec][word] += 1;
            }
            ec += 1;
        }
    }
    return (allWords, articleWords, articleTitles);
}

方法的返回的元组由三部分组成：
allWords统计了每个单词在所有文章出现的次数，该变量用来指示哪些单词应被作为特征的一部分。
articleWords是单词在每篇文章出现的次数。
articleTitles是文章标题的列表。

生成矩阵

接着我们利用上面的字典来生成矩阵。每行矩阵表示了单词在一篇文章出现的次数，而矩阵的列就是对应的单词。
为了使矩阵不至于过大，去掉只存在于个别文章中的单词，及在大部分文章都存在的单词。
实现以上所描述的功能比较简单，在NewsFeatures加入MakeMatrix方法：

public (List<List<double>>, List<string>) MakeMatrix(Dictionary<string, int> allW,
    Dictionary<int, Dictionary<string, int>> articleW)
{
    var wordVec = new List<string>();

    // 只考虑普通的但又不是非常普通的单词
    foreach (var kvp in allW)
    {
        var w = kvp.Key;
        var c = kvp.Value;
        if (c > 3 && c < articleW.Count * 0.6)
            wordVec.Add(w);
    }

    // 构造单词矩阵
    var ll = new List<List<double>>();
    foreach (var f in articleW.Values)
    {
        var r = wordVec.Select(word =>
             {
                 if (f.ContainsKey(word)) return f[word];
                 return 0d;
             })
            .ToList();
        ll.Add(r);
    }
    return (ll, wordVec);
}

这一部分生成的矩阵使用一个嵌套的List来表示。与矩阵同时返回的是单词的列表（即矩阵每一列对应的单词）。
这时就可以测试下解析RSS并生成矩阵（后续所有测试目的的代码都在NewsFeaturesTest类中）：

public class NewsFeaturesTest
{
    private readonly ITestOutputHelper _output;

    public NewsFeaturesTest(ITestOutputHelper output)
    {
        _output = output;
    }

    private void TestOutput(object obj)
    {
        _output.WriteLine(obj.ToString());
    }

    [Fact]
    public void TestMakeMatrix()
    {
        var newsFeatures = new NewsFeatures();
        (var allW,var artW,var artT) = newsFeatures.GetArticleWords();
        (var wordMatrix,var wordVec) = newsFeatures.MakeMatrix(allW, artW);
        // 单词向量的前10个词
        TestOutput(JsonConvert.SerializeObject(wordVec.Take(10)));
        // 第二篇文章的标题
        TestOutput(artT[1]);
        // 第二篇文章的在单词矩阵中对应的数据行的前10个值
        TestOutput(JsonConvert.SerializeObject(wordMatrix[1].Take(10)));
    }
}

非负矩阵因式分解

这是本文最关键的部分。非负矩阵因式分解（NMF）是从数据中提取重要特征的技术。
本文不讲解这个技术的来龙去脉，只从过程分析下将NMF应用于特征提取的过程。

最终我们要得到的结果存在于两个矩阵 - 分别是特征矩阵和权重矩阵。这两个矩阵是由我们之前生成的单词矩阵分解而来，即两个矩阵相乘会得到单词矩阵。

权重矩阵：此矩阵的每一行对应一篇文章，每一列对应一个特征，行列交叉的数字就表示特征与对于文章相关的程度。
特征矩阵：此矩阵中每一行对应一个特征，每一列表示一个单词，行列交叉的数字就表示单词对于对应特征的重要程度。

按照矩阵乘法的规则，将权重矩阵与特征矩阵相乘就可以得到原来的单词矩阵。
这面分解的特征和单词权重都是非负值，这也是非负矩阵因式分解名称来源。

这个要寻找的特征的数量，需要在进行因式分解的时候传入，如果要给每篇文章都找到一个特征，那就把特征的数量指定为文章的数量。如果想要给文章进行一个分类，即找到有共同特征的文章，则把这个特征数量指定为小于文章的数量。
特别注意，现实中很难分解出能完美还原文章的矩阵的两个子矩阵，都是尽量去拟合。从后面的测试可以看出。

关于矩阵分解过程的计算，原书使用Python著名的库NumPy来完成，博主在编写测试时也尝试借助pythonnet在C#使用Python中的NumPy，不过由于类型转换过程的种种问题放弃。后来经过一番寻找发现了.NET平台原生的数学计算库Math.Net。

关于Math.Net的使用可以参见园内大神asxinyu的这个系列文章。本文的测试代码的编写也用到了这些文章中的一部分示例代码。

首先使用Nuget安装Math.NET

Install-Package MathNet.Numerics -Version 3.19.0

然后我们编写一些测试代码，了解Math.NET的使用。

[Fact]
public void TestMathNet()
{
    // 测试矩阵乘法
    var m1 = DenseMatrix.OfArray(new[,] { { 1d, 2d, 3d }, { 4d, 5d, 6d } });
    TestOutput(m1);
    var m2 = DenseMatrix.OfArray(new[,] { { 1d, 2d }, { 3d, 4d }, { 5d, 6d } });
    TestOutput(m2);
    TestOutput((m1 * m2).ToString());

    //测试生成随机数组
    var mb = Matrix<double>.Build;
    var rnd = new Random();
    var mr = mb.Dense(2, 3, (i, j) => rnd.NextDouble());
    TestOutput(mr);

    //测试矩阵转一维数组并由数组重构矩阵
    var arr = mr.AsColumnMajorArray();
    TestOutput(JsonConvert.SerializeObject(arr));

    var mback = mb.DenseOfColumnMajor(2, 3, arr);
    TestOutput(mback);

    //测试由行列表生成矩阵，并转换
    var rowList = new List<List<double>>()
    {
        new List<double>(){1, 2,3,4 },
        new List<double>(){5,6,7,8 },
        new List<double>(){9,10,11,12}
    };
    var matrixFromRow = mb.DenseOfRows(rowList);
    var matrixFromRowToColumnArr = matrixFromRow.AsColumnMajorArray();
    TestOutput(JsonConvert.SerializeObject(matrixFromRowToColumnArr));
}

这个方法中还测试了将矩阵转为一维数组，并由一维数组还原矩阵，具体用途后文可以看到。

NMF算法

具体的NMF算法的证明也很复杂，这里也只是描述下过程，并直接给出实现。

首先需要一个成本函数来衡量分解的效果，方法就是比较分解得到的子矩阵相乘后的结果与原矩阵的差异。而比较两个矩阵的方法是计算矩阵中每个值的差的平方的和。新建一个Nmf类，在其中添加DifCost方法：

public class Nmf
{
    private readonly Action<object> _outputWriter;

    public Nmf(Action<object> outputAction)
    {
        _outputWriter = outputAction;
    }

    private double DifCost(DenseMatrix a, Matrix<double> b)
    {
        var dif = 0d;
        // 遍历矩阵中的每一行和每一列
        for (int i = 0; i < a.RowCount; i++)
        {
            for (int j = 0; j < a.ColumnCount; j++)
            {
                //将差值相加
                dif += Math.Pow(a[i, j] - b[i, j], 2);
            }
        }
        return dif;
    }
}

接着借助乘法更新法则(multiplicative update rules)来逐步更新矩阵，使上面的成本函数的计算值逐步降低。
乘法更新法则也是谜一般的存在，其方法如下：

法则的计算产生4个新的更新矩阵。

• hn 经过转置后的权重矩阵与单词矩阵相乘得到的矩阵
• hd 经过转置后的权重矩阵与原权重矩阵相乘，再与特征矩阵相乘得到的矩阵
• wn 单词矩阵与经过转置后的特征矩阵相乘得到的矩阵
• wd 权重矩阵与特征矩阵相乘，再与转置后的特征矩阵相乘得到的矩阵

为了更新特征矩阵，首先将上述所有矩阵转为数组，将特征矩阵的每一个值与hn中的对应值相乘，并除以hd中的对应值。权重矩阵的更新方法类似，见算法实现。
在Nmf类中添加Factorize方法来进行分解计算。
其中，w表示权重矩阵，h为特征矩阵。
方法接受的参数中pc，表示要寻找的特征的数量，即权重矩阵w的列数，也即特征矩阵h的行数。

特征数量的选取，有的业务场景下有明确的要求，而有时候据需要更具经验，或反复的实验来选取合适的数目。

public (Matrix<double>, Matrix<double>) Factorize(DenseMatrix v, int pc = 10, int iter = 50)
{
    var ic = v.RowCount;
    var fc = v.ColumnCount;

    //以随机值初始化权重矩阵和特征矩阵
    var mb = Matrix<double>.Build;
    var rndGen = new Random(DateTime.Now.Second);
    var w = mb.Dense(ic,pc, (i, j) => rndGen.NextDouble());
    var h = mb.Dense(pc,fc, (i, j)=> rndGen.NextDouble());

    // 最多执行iter次操作
    for (int i = 0; i < iter; i++)
    {
        var wh = w * h;
        // 计算当前查值
        var cost = DifCost(v, wh);
        if (i % 10 == 0)
            _outputWriter(cost);
        //如果矩阵已分解彻底，终止循环
        if (cost == 0)
            break;
        //更新特征矩阵
        var hn = w.Transpose() * v;
        var hd = w.Transpose() * w * h;
        h = DenseMatrix.OfColumnMajor(pc, fc,
            h.AsColumnMajorArray().Multiple(hn.AsColumnMajorArray()).Divide(hd.AsColumnMajorArray()));

        //更新权重矩阵
        var wn = v * h.Transpose();
        var wd = w * h * h.Transpose();
        w = DenseMatrix.OfColumnMajor(ic, pc,
            w.AsColumnMajorArray().Multiple(wn.AsColumnMajorArray()).Divide(wd.AsColumnMajorArray()));
    }
    return (w, h);
}

Math.Net不支持原书中所需要的将矩阵作为数组并进行相乘相除的操作（也可能是楼主没找到），所以实现了两个扩展方法用于数组的乘除操作。为了避免除0还做了一点小小的处理。

public static class DoubleArrayExt
{
    public static double[] Multiple(this double[] left, double[] right)
    {
        return left.Select((l, i) => l * right[i]).ToArray();
    }

    public static double[] Divide(this double[] left, double[] right)
    {
        return left.Select((l, i) => l / (right[i]==0?double.MinValue:right[i])).ToArray();
    }
}

Nmf算法的实现就是这些，可以尝试来进行下矩阵分解：

[Fact]
public void TestFactorize()
{
    var m1 = DenseMatrix.OfArray(new[,] { { 1d, 2d, 3d }, { 4d, 5d, 6d } });
    var m2 = DenseMatrix.OfArray(new[,] { { 1d, 2d }, { 3d, 4d }, { 5d, 6d } });
    var m12 = m1 * m2;
    TestOutput(m12);

    var nmf = new Nmf(o => TestOutput(o.ToString()));
    (var w, var h) = nmf.Factorize(m12, 3, 100);
    TestOutput(w);
    TestOutput(h);
    TestOutput(w * h);
}

注意每次分解所得的矩阵都不同，但分解的矩阵相乘都能得到原矩阵。

单词矩阵的分解及特征的提取

有了单词矩阵和分解算法，下面就可以尝试将单词矩阵分解：

[Fact]
public void TestArticleFactorize()
{
    var newsFeatures = new NewsFeatures();
    (var allW, var artW, var artT) = newsFeatures.GetArticleWords();
    (var wordMatrix, var wordVec) = newsFeatures.MakeMatrix(allW, artW);
    var nmf = new Nmf(o => TestOutput(o.ToString()));
    var matrix = DenseMatrix.OfRows(wordMatrix);
    (var weights, var feat) = nmf.Factorize(matrix, 20, 50);
    TestOutput(weights);
    TestOutput(feat);
}

结果中，weights就是保存特征应用于文章的权重的矩阵，而feat是每篇文章的特征（单词对特征贡献度）矩阵。
通过计算过程输出的成本函数计算结果可以看出，分解出的矩阵（相乘后）与原矩阵拟合度不是太高。

下面通过一种更直观的方法来展示权重矩阵和特征矩阵的所体现的结果。
首先我们的特征是以单词应用于特征的权重来表示，我们由特征矩阵提取每个对每个特征贡献度最高的前6个单词来表示一个特征。
然后我们看看与一个特征最相关的前3篇文章，在NewsFeature类中添加ShowFeatures方法：

public (List<List<ArticleAndFeature>>, List<string>) ShowFeatures(Matrix<double> w, Matrix<double> h,
    List<string> titles, List<string> wordvec, string @out = "features.txt")
{
    using (var fs = File.Create(@out))
    {
        var sw = new StreamWriter(fs);
        var pc = h.RowCount;
        var wc = h.ColumnCount;
        var topPatterns = Enumerable.Repeat(new List<ArticleAndFeature>(), titles.Count).ToList();
        var patternNames = new List<string>();

        //遍历所有特征
        for (int i = 0; i < pc; i++)
        {
            var slist = new List<ValueTuple<double, string>>();
            //构造包含单词及其权重数据的列表
            for (int j = 0; j < wc; j++)
            {
                slist.Add((h[i, j], wordvec[j]));
            }
            //按单词对特征贡献度值倒叙排序
            slist.Sort((x,y)=>y.Item1.CompareTo(x.Item1));
            //打印开始的6个元素
            var n = string.Join(",", slist.Take(6));
            sw.WriteLine($"[{n}]");
            patternNames.Add(n);

            //构造针对该特征的文章列表
            var flist = new List<ValueTuple<double, string>>();
            for (int j = 0; j < titles.Count; j++)
            {
                //加入文章及权重数据
                flist.Add((w[j, i], titles[j]));
                topPatterns[j].Add(new ArticleAndFeature(w[j, i], i, titles[j]));
            }
            // 按特征对于文章的匹配程度有高到低排列
            flist.Sort((x,y)=>y.Item1.CompareTo(x.Item1));
            //显示前3篇文章
            foreach (var kvp in flist.Take(3))
            {
                sw.WriteLine($"({kvp.Item1},{kvp.Item2})");
                sw.WriteLine();
            }
        }
        sw.Close();
        // 返回模式名称，后面要用
        return (topPatterns, patternNames);
    }
}

上面的方法中用了一个保存文章和特征相关度的子类ArticleAndFeature，也将其加入NewsFeatures中

public class ArticleAndFeature : IComparable<ArticleAndFeature>
{
    public ArticleAndFeature(double weight, int featureIdx, string articleTitle)
    {
        Weight = weight;
        FeatureIdx = featureIdx;
        ArticleTitle = articleTitle;
    }

    public double Weight { get; set; }

    public int FeatureIdx { get; set; }

    public string ArticleTitle { get; set; }

    public int CompareTo(ArticleAndFeature right)
    {
        return right.Weight.CompareTo(Weight);//由大到小排序
    }
}

由于输出的结果会非常长，示例代码中将结果保存到文件里。
然后可以用下面的代码来测试下看看与特征密切度最高的3篇文章。

[Fact]
public void TestShowFeatures()
{
    var newsFeatures = new NewsFeatures();
    (var allW, var artW, var artT) = newsFeatures.GetArticleWords();
    (var wordMatrix, var wordVec) = newsFeatures.MakeMatrix(allW, artW);
    var nmf = new Nmf(o => TestOutput(o.ToString()));
    var matrix = DenseMatrix.OfRows(wordMatrix);
    (var weights, var feat) = nmf.Factorize(matrix, 20, 50);
    newsFeatures.ShowFeatures(weights, feat, artT, wordVec);
}

之前我们列出与一个特征最相关的三篇文章，反过来也可以查看与一篇文章最相关的三个特征，而且从特征的相关度数值看出文章到底只有一个主题(特征)还是与列出的主题都有密切的关系。方法ShowArticles也加入NewsFeatures中。

public void ShowArticles(List<string> titles,
    List<List<ArticleAndFeature>> topPatterns,
    List<string> patternNames, string @out = "articles.txt")
{
    using (var fs = File.Create(@out))
    {
        var sw = new StreamWriter(fs);

        // 遍历所有文章
        for (int j = 0; j < titles.Count; j++)
        {
            sw.WriteLine(titles[j]);

            // 针对当前文章，获得排位最靠前(倒序下)的几个特征
            topPatterns[j].Sort();

            // 打印前3个模式
            for (int i = 0; i < 3; i++)
            {
                sw.WriteLine($@"{topPatterns[j][i].Weight}
                    {JsonConvert.SerializeObject(patternNames[topPatterns[j][i].FeatureIdx])}");
            }
            sw.WriteLine();
        }
        sw.Close();
    }
}

可以用下面的代码来测试，并查看输出到文件的结果：

[Fact]
public void TestShowArticles()
{
    var newsFeatures = new NewsFeatures();
    (var allW, var artW, var artT) = newsFeatures.GetArticleWords();
    (var wordMatrix, var wordVec) = newsFeatures.MakeMatrix(allW, artW);
    var nmf = new Nmf(o => TestOutput(o.ToString()));
    var matrix = DenseMatrix.OfRows(wordMatrix);
    (var weights, var feat) = nmf.Factorize(matrix, 20, 50);
    (var topPatterns, var patternNames) = newsFeatures.ShowFeatures(weights, feat, artT, wordVec);
    newsFeatures.ShowArticles(artT, topPatterns, patternNames);
}

这一节的主要内容就到此结束，我们了解到了如果使用NMF方法进行特征提取。

附

楼主曾尝试在测试代码中使用pythonnet作为桥梁来使用Numpy进行矩阵计算，虽然未果，但还用把pythonnet的使用记录下来给有需要的园友参考。

关于在C#中使用NumPy可以使用这个库。这个库不像IronPython那样在CLR上完整的实现了一遍Python，而是通过加载一个CPython来实现Python库及代码的执行。
虽然nuget上有发布的pythonnet包，但只支持32位的Python环境，而最重要的是这个包不知道由于什么问题不支持numpy的数组和矩阵计算。
所以我们自行下载源代码来编译，写本文时楼主下载的pythonnet的提交版本为550a027。
由于楼主电脑安装了Python3.6.1，使用VS打开解决方案，将生成配置改为DebugWinPY3或ReleaseWinPY3之一，如图：

这两个配置下的条件编译符号是我们所需要的：

把生成的Python.Runtime.dll加入到我们的测试项目中待用。

最新的3.6.1版的64位python安装版，官网下载地址。使用安装包的好处是其可以注册到系统环境中。如果没有这个注册过程pythonnet也无法正确加载python运行环境。

在安装Python后要重启一下计算机，pythonnet才能正常的将python运行环境加载到CLR中，否则会报python36未找到等类似错误。

官方的python安装包是没有集成numpy模块的，不过好在目前的python安装包都集成了pip。
在Windows命令行下执行pip install numpy安装numpy，安装后进入python交互环境，使用help(modules)确认模块列表中已经有了numpy模块。

另外pythonnet库需要.NET Framework4.0以上版本，测试项目使用了.NET Framework4.7符合要求，如果是调用64位平台生成的Python.Runtime项目，调用者（即PCI项目）也必须是64位。

环境准备好后，运行下面的测试代码看看是否可以正常的进行矩阵乘法计算。如果计算说明一切正常。可以继续后面的工作。

[Fact]
public void TestNumpy()
{
    using (Py.GIL())
    {
        Numpy.Initialize();
        var scope = Py.CreateScope();
        var np = scope.Import("numpy", "np");


        dynamic m1 = np.matrix(new List<List<float>> { new List<float>() { 1, 2, 3 }, new List<float>() { 4, 5, 6 } });
        TestOutput(m1);
        dynamic m2 = np.matrix(new List<List<float>> { new List<float>() { 1, 2,  }, new List<float>() { 3,4},new List<float>() { 5, 6 } });
        TestOutput(m2);

        TestOutput((m1*m2).ToString());
    }
}

本系列示例代码下载。