《集體智慧編程》讀書筆記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());
}
}
本系列示例代碼下載。
浙公網安備 33010602011771號