Permalink
Switch branches/tags
Nothing to show
Find file Copy path
Fetching contributors…
Cannot retrieve contributors at this time
195 lines (162 sloc) 7.4 KB
(*
Implementation of the Non-Negative Matrix Factorization algorithm in Mathematica
Copyright (C) 2013 Anton Antonov
This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 3 of the License, or
(at your option) any later version.
This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
GNU General Public License for more details.
You should have received a copy of the GNU General Public License
along with this program. If not, see <http://www.gnu.org/licenses/>.
Written by Anton Antonov,
antononcube @ gmail . com,
Windermere, Florida, USA.
*)
(*
Mathematica is (C) Copyright 1988-2013 Wolfram Research, Inc.
Protected by copyright law and international treaties.
Unauthorized reproduction or distribution subject to severe civil
and criminal penalties.
Mathematica is a registered trademark of Wolfram Research, Inc.
*)
(* Version 1.0 *)
(* This package contains definitions for the application of Non-Negative Matrix Factorization (NNMF). *)
(*
The implementation follows the description of the hybrid algorithm GD-CLS (Gradient Descent with Constrained Least Squares) in the article:
Shahnaz, F., Berry, M., Pauca, V., Plemmons, R., 2006.
Document clustering using nonnegative matrix factorization. Information Processing & Management 42 (2), 373-386.
In order to use NearestWords a nearest function has to be created over the column indices of the topic matrix.
For example:
{W, H} = NormalizeMatrixProduct[W, H];
HNF = Nearest[Range[Dimensions[H][[2]]], DistanceFunction -> (Norm[H[[All, #1]] - H[[All, #2]]] &)]
NearestWords[HNF, "agent", termsOfH, stemmingRules, 15]
*)
BeginPackage["NonNegativeMatrixFactorization`"]
GDCLS::usage = "GDCLS[V_?MatrixQ,k_Integer,opts] returns the pair of matrices {W,H} such that V = W H and \
the number of the columns of W and the number of rows of H are k. The method used is called Gradient Descent with Constrained Least Squares."
GDCLSGlobal::usage = "GDCLSGlobal[V_?MatrixQ,W_?MatrixQ,H_?MatrixQ,opts] continues the GDCLS iterations over the matrices W and H \
in the execution context and returns {W,H} as a result."
NormalizeMatrixProduct::usage = "NormalizeMatrixProduct[W_?MatrixQ,H_?MatrixQ] returns a pair of matrices {W1,H1} \
such that W1 H1 = W H and the norms of the columns of W1 are 1."
LeftNormalizeMatrixProduct::usage = "Same as NormalizeMatrixProduct."
RightNormalizeMatrixProduct::usage = "RightNormalizeMatrixProduct[W_?MatrixQ,H_?MatrixQ] returns a pair of matrices {W1,H1} \
such that W1 H1 = W H and the norms of the rows of H1 are 1."
BasisVectorInterpretation::usage = "BasisVectorInterpretation[vec_?VectorQ,n_Integer,interpretationItems_List] \
takes the n largest coordinates of vec, finds the corresponding elements in interpretationItems, and returns a list of coordinate-item pairs."
NearestWords::usage = "NearestWords[HNF, word, terms, stemmingRules, n] calculates a statistical thesaurus entry \
for a specified nearest function over the columns of a matrix of topics and a word."
Begin["`Private`"]
Clear[GDCLS]
Options[GDCLS] = {"MaxSteps" -> 200, "NonNegative" -> True, "Epsilon" -> 10^-9., "RegularizationParameter" -> 0.01, PrecisionGoal -> Automatic, "PrintProfilingInfo" -> False};
GDCLS[V_?MatrixQ, k_?IntegerQ, opts:OptionsPattern[]] :=
Block[{t, fls, A, W, H, T, m, n, b, diffNorm, normV, nSteps = 0,
nonnegQ = OptionValue[GDCLS,"NonNegative"],
maxSteps = OptionValue[GDCLS,"MaxSteps"],
eps = OptionValue[GDCLS,"Epsilon"],
lbd = OptionValue[GDCLS,"RegularizationParameter"],
pgoal = OptionValue[GDCLS,PrecisionGoal],
PRINT = If[TrueQ[OptionValue[GDCLS,"PrintProfilingInfo"]], Print, None]},
{m, n} = Dimensions[V];
W = RandomReal[{0, 1}, {m, k}];
H = ConstantArray[0, {k, n}];
normV = Norm[V, "Frobenius"]; diffNorm = 10 normV;
If[ pgoal === Automatic, pgoal = 4 ];
While[nSteps < maxSteps && TrueQ[! NumberQ[pgoal] || NumberQ[pgoal] && (normV > 0) && diffNorm/normV > 10^(-pgoal)],
nSteps++;
t =
Timing[
A = Transpose[W].W + lbd*IdentityMatrix[k];
T = Transpose[W];
fls = LinearSolve[A];
H = Table[(b = T.V[[All, i]]; fls[b]), {i, 1, n}];
H = SparseArray[Transpose[H]];
If[nonnegQ,
H = Clip[H, {0, Max[H]}]
];
W = W*(V.Transpose[H])/(W.(H.Transpose[H]) + eps);
];
If[NumberQ[pgoal],
diffNorm = Norm[V - W.H, "Frobenius"];
If[nSteps < 100 || Mod[nSteps, 100] == 0, PRINT["step:", nSteps, ", iteration time:", t, " relative error:", diffNorm/normV]],
If[nSteps < 100 || Mod[nSteps, 100] == 0, PRINT["step:", nSteps, ", iteration time:", t]]
];
];
{W, H}
];
Clear[GDCLSGlobal]
Options[GDCLSGlobal] = Options[GDCLS];
SetAttributes[GDCLSGlobal,HoldAll]
GDCLSGlobal[V_, W_, H_, opts:OptionsPattern[]] :=
Block[{t, fls, A, k, T, m, n, b, diffNorm, normV, nSteps = 0,
nonnegQ = OptionValue[GDCLSGlobal,"NonNegative"],
maxSteps = OptionValue[GDCLSGlobal,"MaxSteps"],
eps = OptionValue[GDCLSGlobal,"Epsilon"],
lbd = OptionValue[GDCLSGlobal,"RegularizationParameter"],
pgoal = OptionValue[GDCLSGlobal,PrecisionGoal],
PRINT = If[TrueQ[OptionValue[GDCLSGlobal,"PrintProfilingInfo"]], Print, None]},
{m, n} = Dimensions[V];
k = Dimensions[H][[1]];
normV = Norm[V, "Frobenius"]; diffNorm = 10 normV;
While[nSteps < maxSteps && TrueQ[! NumberQ[pgoal] || NumberQ[pgoal] && (normV > 0) && diffNorm/normV > 10^(-pgoal)],
nSteps++;
t =
Timing[
A = Transpose[W].W + lbd*IdentityMatrix[k];
T = Transpose[W];
fls = LinearSolve[A];
H = Table[(b = T.V[[All, i]]; fls[b]), {i, 1, n}];
H = SparseArray[Transpose[H]];
If[nonnegQ,
H = Clip[H, {0, Max[H]}]
];
W = W*(V.Transpose[H])/(W.(H.Transpose[H]) + eps);
];
If[NumberQ[pgoal],
diffNorm = Norm[V - W.H, "Frobenius"];
If[nSteps < 100 || Mod[nSteps, 100] == 0, PRINT[nSteps, " ", t, " relative error=", diffNorm/normV]],
If[nSteps < 100 || Mod[nSteps, 100] == 0, PRINT[nSteps, " ", t]]
];
];
{W, H}
] /; MatrixQ[W] && MatrixQ[H] && Dimensions[W][[2]] == Dimensions[H][[1]];
(* ::Subsection:: *)
(*Normalize matrices*)
Clear[NormalizeMatrixProduct]
NormalizeMatrixProduct[W_?MatrixQ,H_?MatrixQ]:=
Block[{d,S,SI},
d=Table[Norm[W[[All,i]]],{i,Length[W[[1]]]}];
S=DiagonalMatrix[d];
SI=DiagonalMatrix[Map[If[#!=0,1/#,0]&,d]];
{W.(SI),S.H}
];
LeftNormalizeMatrixProduct = NormalizeMatrixProduct;
Clear[RightNormalizeMatrixProduct]
RightNormalizeMatrixProduct[W_?MatrixQ,H_?MatrixQ]:=
Block[{d,S,SI},
d=Table[Norm[H[[i]]],{i,Length[H]}];
S=DiagonalMatrix[d];
SI=DiagonalMatrix[1/d];
{W.S,SI.H}
];
Clear[BasisVectorInterpretation]
BasisVectorInterpretation[vec_,n_Integer,terms_]:=
Block[{t},
t=Reverse@Ordering[vec,-n];
Transpose[{vec[[t]],terms[[t]]}]
];
Clear[NearestWords];
NearestWords[HNF_NearestFunction, word_String, terms : {_String ..},
stemmingRules_, n_Integer: 20] :=
Block[{sword, tpos, inds},
sword = word /. stemmingRules;
tpos = Position[terms, sword];
If[Length[tpos] == 0, {},
inds = HNF[Flatten[tpos][[1]], n];
terms[[inds]]
]
];
End[]
EndPackage[]