d-spca.html

buildbot Fri, 03 Feb 2017 14:26:25 -0800

Author: buildbot
Date: Fri Feb  3 22:26:05 2017
New Revision: 1006162

Log:
Staging update by buildbot for mahout


Modified:
    websites/staging/mahout/trunk/content/   (props changed)
    websites/staging/mahout/trunk/content/users/algorithms/d-spca.html

Propchange: websites/staging/mahout/trunk/content/
------------------------------------------------------------------------------
--- cms:source-revision (original)
+++ cms:source-revision Fri Feb  3 22:26:05 2017
@@ -1 +1 @@
-1781487
+1781612

Modified: websites/staging/mahout/trunk/content/users/algorithms/d-spca.html
==============================================================================
--- websites/staging/mahout/trunk/content/users/algorithms/d-spca.html 
(original)
+++ websites/staging/mahout/trunk/content/users/algorithms/d-spca.html Fri Feb  
3 22:26:05 2017
@@ -281,43 +281,163 @@
 h2:hover > .headerlink, h3:hover > .headerlink, h1:hover > .headerlink, 
h6:hover > .headerlink, h4:hover > .headerlink, h5:hover > .headerlink, 
dt:hover > .elementid-permalink { visibility: visible }</style>
 <h1 id="distributed-stochastic-pca">Distributed Stochastic PCA<a 
class="headerlink" href="#distributed-stochastic-pca" title="Permanent 
link">&para;</a></h1>
 <h2 id="intro">Intro<a class="headerlink" href="#intro" title="Permanent 
link">&para;</a></h2>
-<p>Mahout has a distributed implementation of Stochastic PCA </p>
-<h2 id="motivation">Motivation<a class="headerlink" href="#motivation" 
title="Permanent link">&para;</a></h2>
-<p>Stochastic SVD method in Mahout produces reduced-rank Singular Value 
Decomposition output in its strict mathematical definition: 
<code>\(\mathbf{A}\approx\mathbf{UÎ£V}\)</code>, i.e. it creates outputs for 
matrices <code>\(\mathbf{U},\mathbf{V}, and \mathbf{Î£}\)</code>, each of which 
may be requested individually. The desired rank of decomposition, henceforth 
denoted as <em>k</em><code>\(\in\mathbb{N}_1\)</code>, is a parameter of the 
algorithm. The singular values inside diagonal matrix <code>\(\Sigma\)</code> 
satisfyÏi+1â¤Ïiâiâ[1,kâ1], i.e. sorted from biggest tosmallest. Cases 
of rank deficiency rank(A)&lt; karehandled by producing 0s in singular value 
positionsonce deficiency takes place.</p>
+<p>Mahout has a distributed implementation of Stochastic PCA[1].</p>
+<h2 id="algorithm">Algorithm<a class="headerlink" href="#algorithm" 
title="Permanent link">&para;</a></h2>
+<p>Given an <em>m</em> <code>\(\times\)</code> <em>n</em> matrix 
<code>\(\mathbf{A}\)</code>, a target rank <em>k</em>, and an oversampling 
parameter <em>p</em>, this procedure computes a <em>k</em>-rank PCA by finding 
the unknowns in <code>\(\mathbf{A}â1\mu \ge U\SigmaV\)</code>:
+(1) Create seed for random <em>n</em> <code>\(\times\)</code> <em>(k+p)</em> 
matrix <code>\(\Omega\)</code>.
+(2) <code>\(s_\Omega \leftarrow \Omega^\top \mu\)</code>.
+(3) <code>\(\mathbf{Y_0} \leftarrow A\Omega â 1(s_\Omega)^\top, Y \in 
\mathbb{R}^(m\times(k+p))\)</code>.
+(4) Column-orthonormalize <code>\(\mathbf{Y_0} \rightarrow \mathbf{Q}\)</code> 
by computing thin decomposition <code>\(\mathbf{Y_0} = \mathbf{QR}\)</code>. 
Also, <code>\(\mathbf{Q}\in\mathbb{R}^(m\times(k+p)), 
\mathbf{R}\in\mathbb{R}^((k+p)\times(k+p))\)</code>.
+(5) <code>\(s_Q \leftarrow Q^\top 1\)</code>.
+(6) <code>\(\mathbf{B_0} \leftarrow Q^\top A: B \in \mathbb{R}^((k+p)\times 
n)\)</code>.
+(7) <code>\(s_B \leftarrow (B_0)^\top \mu\)</code>.
+(8) For <em>i</em> in 1..<em>q</em> repeat (power iterations):
+  (a) For <em>j</em> in 1..<em>n</em> <code>\(apply(B_(iâ1))_(âj) 
\leftarrow (B_(iâ1))_(âj)â\mu_j s_Q\)</code>.
+  (b) <code>\(\mathbf{Y_i) \leftarrow 
\mathbf{(AB_(iâ1)^\top)â1(s_Bâ\mu^\top \mu s_Q^\top)}\)</code>.
+  (c) Column-orthonormalize <code>\(\mathbf{Y_i} \rightarrow 
\mathbf{Q}\)</code> by computing thin decomposition <code>\(\mathbf{Y_i = 
QR}\)</code>.
+  (d) <code>\(\mathbf{s_Q \leftarrow Q^\top 1}\)</code>.
+  (e) <code>\(\mathbf{B_i \leftarrow Q^\top A}\)</code>.
+  (f) <code>\(\mathbf{s_B \leftarrow (B_i)^\top \mu}\)</code>.
+(9) Let <code>\(\mathbf{C \triangleq s_Q (s_B)^\top}\)</code>. 
<code>\(\mathbf{M \leftarrow B_q (B_q)^\top â C â C^\top + \mu^\top \mu s_Q 
(s_Q)^\top}\)</code>.
+(10) Compute an eigensolution of the small symmetric <code>\(\mathbf{M = 
\hat{U} \Lambda \hat{U}^\top: M \in \mathbb{R}^((k+p)\times(k+p))}\)</code>.
+(11) The singular values <code>\(\Sigma = \Lambda^(\circ 0.5)\)</code>, or, in 
other words, <code>\(\mathbf{\sigma_i= \sqrt{\lambda_i}}\)</code>.
+(12) If needed, compute <code>\(\mathbf{U = Q\hat{U}}\)</code>.
+(13) If needed, compute <code>\(\mathbf{V = B^\top \hat{U} 
\Sigma^(â1)}\)</code>. Another way is <code>\(\mathbf{V = A^\top 
U\Sigma^(â1)}\)1.
+(14) If needed, items converted to the PCA space can be computed 
as</code>(\mathbf{U\Sigma})`.</p>
 <h2 id="implementation">Implementation<a class="headerlink" 
href="#implementation" title="Permanent link">&para;</a></h2>
-<p>Mahout <code>dqrThin(...)</code> is implemented in the mahout 
<code>math-scala</code> algebraic optimizer which translates Mahout's R-like 
linear algebra operators into a physical plan for both Spark and H2O 
distributed engines.</p>
-<div class="codehilite"><pre><span class="n">def</span> <span 
class="n">dqrThin</span><span class="p">[</span><span class="n">K</span><span 
class="p">:</span> <span class="n">ClassTag</span><span 
class="p">](</span><span class="n">A</span><span class="p">:</span> <span 
class="n">DrmLike</span><span class="p">[</span><span class="n">K</span><span 
class="p">],</span> <span class="n">checkRankDeficiency</span><span 
class="p">:</span> <span class="n">Boolean</span> <span class="p">=</span> 
<span class="n">true</span><span class="p">):</span> <span 
class="p">(</span><span class="n">DrmLike</span><span class="p">[</span><span 
class="n">K</span><span class="p">],</span> <span class="n">Matrix</span><span 
class="p">)</span> <span class="p">=</span> <span class="p">{</span>        
-    <span class="k">if</span> <span class="p">(</span><span 
class="n">drmA</span><span class="p">.</span><span class="n">ncol</span> <span 
class="o">&gt;</span> 5000<span class="p">)</span>
-        <span class="nb">log</span><span class="p">.</span><span 
class="n">warn</span><span class="p">(</span>&quot;<span class="n">A</span> 
<span class="n">is</span> <span class="n">too</span> <span 
class="n">fat</span><span class="p">.</span> <span class="n">A</span><span 
class="o">&#39;</span><span class="n">A</span> <span class="n">must</span> 
<span class="n">fit</span> <span class="n">in</span> <span 
class="n">memory</span> <span class="n">and</span> <span 
class="n">easily</span> <span class="n">broadcasted</span><span 
class="p">.</span>&quot;<span class="p">)</span>
-    <span class="n">implicit</span> <span class="n">val</span> <span 
class="n">ctx</span> <span class="p">=</span> <span class="n">drmA</span><span 
class="p">.</span><span class="n">context</span>
-    <span class="n">val</span> <span class="n">AtA</span> <span 
class="p">=</span> <span class="p">(</span><span class="n">drmA</span><span 
class="p">.</span><span class="n">t</span> <span class="c">%*% 
drmA).checkpoint()</span>
-    <span class="n">val</span> <span class="n">inCoreAtA</span> <span 
class="p">=</span> <span class="n">AtA</span><span class="p">.</span><span 
class="n">collect</span>
-    <span class="n">val</span> <span class="n">ch</span> <span 
class="p">=</span> <span class="n">chol</span><span class="p">(</span><span 
class="n">inCoreAtA</span><span class="p">)</span>
-    <span class="n">val</span> <span class="n">inCoreR</span> <span 
class="p">=</span> <span class="p">(</span><span class="n">ch</span><span 
class="p">.</span><span class="n">getL</span> <span 
class="n">cloned</span><span class="p">)</span> <span class="n">t</span>
-    <span class="k">if</span> <span class="p">(</span><span 
class="n">checkRankDeficiency</span> <span class="o">&amp;&amp;</span> !<span 
class="n">ch</span><span class="p">.</span><span 
class="n">isPositiveDefinite</span><span class="p">)</span>
-        <span class="n">throw</span> <span class="n">new</span> <span 
class="n">IllegalArgumentException</span><span class="p">(</span>&quot;<span 
class="n">R</span> <span class="n">is</span> <span class="n">rank</span><span 
class="o">-</span><span class="n">deficient</span><span 
class="p">.</span>&quot;<span class="p">)</span>
-    <span class="n">val</span> <span class="n">bcastAtA</span> <span 
class="p">=</span> <span class="n">sc</span><span class="p">.</span><span 
class="n">broadcast</span><span class="p">(</span><span 
class="n">inCoreAtA</span><span class="p">)</span>
-    <span class="n">val</span> <span class="n">Q</span> <span 
class="p">=</span> <span class="n">A</span><span class="p">.</span><span 
class="n">mapBlock</span><span class="p">()</span> <span class="p">{</span>
-        <span class="k">case</span> <span class="p">(</span><span 
class="n">keys</span><span class="p">,</span> <span class="n">block</span><span 
class="p">)</span> <span class="p">=</span><span class="o">&gt;</span> <span 
class="n">keys</span> <span class="o">-&gt;</span> <span 
class="n">chol</span><span class="p">(</span><span 
class="n">bcastAtA</span><span class="p">).</span><span 
class="n">solveRight</span><span class="p">(</span><span 
class="n">block</span><span class="p">)</span>
-    <span class="p">}</span>
-    <span class="n">Q</span> <span class="o">-&gt;</span> <span 
class="n">inCoreR</span>
-<span class="p">}</span>
+<p>Mahout <code>dspca(...)</code> is implemented in the mahout 
<code>math-scala</code> algebraic optimizer which translates Mahout's R-like 
linear algebra operators into a physical plan for both Spark and H2O 
distributed engines.</p>
+<p>def dspca<a href="drmA: DrmLike[K], k: Int, p: Int = 15, q: Int = 0">K</a>:
+  (DrmLike[K], DrmLike[Int], Vector) = {</p>
+<div class="codehilite"><pre>// Some mapBlock() calls need it
+implicit val ktag =  drmA.keyClassTag
+
+val drmAcp = drmA.checkpoint()
+implicit val ctx = drmAcp.context
+
+val m = drmAcp.nrow
+val n = drmAcp.ncol
+assert(k <span class="err">&lt;</span>= (m min n), &quot;k cannot be greater 
than smaller of m, n.&quot;)
+val pfxed = safeToNonNegInt((m min n) - k min p)
+
+// Actual decomposition rank
+val r = k + pfxed
+
+// Dataset mean
+val mu = drmAcp.colMeans
+
+val mtm = mu dot mu
+
+// We represent Omega by its seed.
+val omegaSeed = RandomUtils.getRandom().nextInt()
+val omega = Matrices.symmetricUniformView(n, r, omegaSeed)
+
+// This done in front in a single-threaded fashion for now. Even though it 
doesn&#39;t require any
+// memory beyond that is required to keep xi around, it still might be 
parallelized to backs
+// for significantly big n and r. TODO
+val s_o = omega.t %*% mu
+
+val bcastS_o = drmBroadcast(s_o)
+val bcastMu = drmBroadcast(mu)
+
+var drmY = drmAcp.mapBlock(ncol = r) {
+  case (keys, blockA) â
+    val s_o:Vector = bcastS_o
+    val blockY = blockA %*% Matrices.symmetricUniformView(n, r, omegaSeed)
+    for (row â 0 until blockY.nrow) blockY(row, ::) -= s_o
+    keys â blockY
+}
+    // Checkpoint Y
+    .checkpoint()
+
+var drmQ = dqrThin(drmY, checkRankDeficiency = false)._1.checkpoint()
+
+var s_q = drmQ.colSums()
+var bcastVarS_q = drmBroadcast(s_q)
+
+// This actually should be optimized as identically partitioned map-side 
A&#39;B since A and Q should
+// still be identically partitioned.
+var drmBt = (drmAcp.t %*% drmQ).checkpoint()
+
+var s_b = (drmBt.t %*% mu).collect(::, 0)
+var bcastVarS_b = drmBroadcast(s_b)
+
+for (i â 0 until q) {
+
+  // These closures don&#39;t seem to live well with outside-scope vars. This 
doesn&#39;t record closure
+  // attributes correctly. So we create additional set of vals for broadcast 
vars to properly
+  // create readonly closure attributes in this very scope.
+  val bcastS_q = bcastVarS_q
+  val bcastMuInner = bcastMu
+
+  // Fix Bt as B&#39; -= xi cross s_q
+  drmBt = drmBt.mapBlock() {
+    case (keys, block) â
+      val s_q: Vector = bcastS_q
+      val mu: Vector = bcastMuInner
+      keys.zipWithIndex.foreach {
+        case (key, idx) â block(idx, ::) -= s_q * mu(key)
+      }
+      keys â block
+  }
+
+  drmY.uncache()
+  drmQ.uncache()
+
+  val bCastSt_b = drmBroadcast(s_b -=: mtm * s_q)
+
+  drmY = (drmAcp %*% drmBt)
+      // Fix Y by subtracting st_b from each row of the AB&#39;
+      .mapBlock() {
+    case (keys, block) â
+      val st_b: Vector = bCastSt_b
+      block := { (_, c, v) â v - st_b(c) }
+      keys â block
+  }
+      // Checkpoint Y
+      .checkpoint()
+
+  drmQ = dqrThin(drmY, checkRankDeficiency = false)._1.checkpoint()
+
+  s_q = drmQ.colSums()
+  bcastVarS_q = drmBroadcast(s_q)
+
+  // This on the other hand should be inner-join-and-map A&#39;B optimization 
since A and Q_i are not
+  // identically partitioned anymore.
+  drmBt = (drmAcp.t %*% drmQ).checkpoint()
+
+  s_b = (drmBt.t %*% mu).collect(::, 0)
+  bcastVarS_b = drmBroadcast(s_b)
+}
+
+val c = s_q cross s_b
+val inCoreBBt = (drmBt.t %*% drmBt).checkpoint(CacheHint.NONE).collect -=:
+    c -=: c.t +=: mtm *=: (s_q cross s_q)
+val (inCoreUHat, d) = eigen(inCoreBBt)
+val s = d.sqrt
+
+// Since neither drmU nor drmV are actually computed until actually used, we 
don&#39;t need the flags
+// instructing compute (or not compute) either of the U,V outputs anymore. 
Neat, isn&#39;t it?
+val drmU = drmQ %*% inCoreUHat
+val drmV = drmBt %*% (inCoreUHat %*% diagv(1 / s))
+
+(drmU(::, 0 until k), drmV(::, 0 until k), s(0 until k))
 </pre></div>
 
 
+<p>}</p>
 <h2 id="usage">Usage<a class="headerlink" href="#usage" title="Permanent 
link">&para;</a></h2>
-<p>The scala <code>dqrThin(...)</code> method can easily be called in any 
Spark or H2O application built with the <code>math-scala</code> library and the 
corresponding <code>Spark</code> or <code>H2O</code> engine module as 
follows:</p>
+<p>The scala <code>dspca(...)</code> method can easily be called in any Spark 
or H2O application built with the <code>math-scala</code> library and the 
corresponding <code>Spark</code> or <code>H2O</code> engine module as 
follows:</p>
 <div class="codehilite"><pre><span class="n">import</span> <span 
class="n">org</span><span class="p">.</span><span class="n">apache</span><span 
class="p">.</span><span class="n">mahout</span><span class="p">.</span><span 
class="n">math</span><span class="p">.</span><span class="n">_</span>
 <span class="n">import</span> <span class="n">decompositions</span><span 
class="p">.</span><span class="n">_</span>
 <span class="n">import</span> <span class="n">drm</span><span 
class="p">.</span><span class="n">_</span>
 
-<span class="n">val</span><span class="p">(</span><span 
class="n">drmQ</span><span class="p">,</span> <span 
class="n">inCoreR</span><span class="p">)</span> <span class="p">=</span> <span 
class="n">dqrThin</span><span class="p">(</span><span 
class="n">drma</span><span class="p">)</span>
+<span class="n">val</span><span class="p">(</span><span 
class="n">drmQ</span><span class="p">,</span> <span 
class="n">inCoreR</span><span class="p">)</span> <span class="p">=</span> <span 
class="n">dqrThin</span><span class="p">(</span><span 
class="n">drma</span><span class="p">,</span> <span class="n">k</span><span 
class="p">,</span> <span class="n">p</span><span class="p">,</span> <span 
class="n">q</span><span class="p">)</span>
 </pre></div>
 
 
+<p>Note the parameter is optional and its default value is zero.</p>
 <h2 id="references">References<a class="headerlink" href="#references" 
title="Permanent link">&para;</a></h2>
-<p>[1]: <a 
href="http://mahout.apache.org/users/sparkbindings/ScalaSparkBindings.pdf";>Mahout
 Scala and Mahout Spark Bindings for Linear Algebra Subroutines</a></p>
-<p>[2]: <a 
href="http://mahout.apache.org/users/sparkbindings/home.html";>Mahout Spark and 
Scala Bindings</a></p>
+<p>[1]: <a 
href="https://www.amazon.com/Apache-Mahout-MapReduce-Dmitriy-Lyubimov/dp/1523775785";>"Apache
 Mahout: Beyond MapReduce; Distributed Algorithm Design", Lyubimov and 
Palumbo</a></p>
    </div>
   </div>     
 </div>

svn commit: r1006162 - in /websites/staging/mahout/trunk/content: ./ users/algorithms/d-spca.html

Reply via email to