Supercomputing Is a Religion But Clouds Are Mundane

Supercomputer is a specialized device. It can do a lot of math, with precision no other platform can do, but it's about all it can have an edge  for.  Have you ever used one as a webserver? Or at least some sort of middle tier realtime application server? There's no practical benefit in doing it that way. Value/cost ratio is not so good. No reason to use demigod's power for your backyard gardening.

By supercomputers I ofcourse only mean hardware-based NUMA architectures; commodity hardware based  simulations with LAN links (such as Beowulf) IMO are really too slow to be seriously considered NUMA type supercomputers.

Clouds are, on the other hand, are as mundane as a toothbrush. Clouds are as versatile as a minivan. And, you can scale them at a fraction of a cost to have as many cpus as supercomputers do. In that sense, paraphrasing the old saying "today's supercomputer is tomorrow's computer", today's supercomputer is today's cloud. The only problem is there's virtually no shared memory space for those CPUs, and so it takes forever to move data from one CPU to another (compared to NUMA). So you can only harness the power of those cpus if you parallelize work in independent splits big enough so that the cloud could be loaded for long enough before the need for sharing split computations arises again. But wouldn't it be cool if that toothbrush could solve as big problems as its supercomputer counterpart with the same amount of cores might? Wouldn't it be nice for a moment to feel a Demigod's power rushing in one's veins while doing that annoying weeding in your backyard?

So, what can we do about highly interconnected problems that rely a lot on shared memory to parallelize a solution?

Fortunately, many interconnected problems in their traditional solutions are connected much more than necessary just because software abstractions such as MPI make shared memory access look cheap and iterative approach is used more than it needs be. Or at least less iterative approach exists. It's just sometimes you need a Demigod's wisdom to see it. But as they say it in the east, the need is the forefather of all inventions.

Follow-up for the "Mean Summarizer..." post

This is a small follow-up for the previous post.

In practice, in MapReduce world, such as Pig UDF functions (which I btw already put in place), when we run over observation data, we often encounter 2 types of problems:


1-- Unordered observations, i.e. cases such as $$t_{n+1}<t_{n}$$.
2-- Parallel processing of disjoint subsets of common superset of observations and combining into one ("combiner"-hadoop, "Algebraic UDF"-pig).

Hence, we need to add those specific cases to our formulas.

1. Updating state with events in the past.
Updated formula will look like:

Case $$t_{n+1}\geq t_{n} $$:
\[\begin{cases}\begin{cases}\pi_{0}=1,\\\pi_{n+1}=e^{-\left(t_{n+1}-t_{n}\right)/\alpha};\end{cases}\\w_{n+1}=1+\pi_{n+1}w_{n};\\s_{n+1}=x_{n+1}+\pi_{n+1}s_{n};\\t_{n+1}=t_{n+1}.\end{cases}\]

Case $$t_{n+1}<t_{n}$$ (updating in-the-past):
\[\begin{cases}\begin{cases}\pi_{0}=1,\\\pi=e^{-\left(t_{n}-t_{n+1}\right)/\alpha},\end{cases}\\w_{n+1}=w_{n}+\pi_{n+1},\\s_{n+1}=s_{n}+\pi_{n+1}x_{n+1},\\t_{n+1}=t_{n}.\end{cases}\]

2.  Combining two summarizers having observed two disjoint sets as subsets of original observation set.
Combining two summarizer states $$S_{1},\, S_{2}$$ having observed two disjoint sets of original observation superset:
Case $$t_{2}\geq t_{1}$$:
\[\begin{cases}t=t_{2};\\s=s_{2}+s_{1}e^{-\left(t_{2}-t_{1}\right)/\alpha};\\w=w_{2}+w_{1}e^{-\left(t_{2}-t_{1}\right)/\alpha}.\end{cases}\]

Case $$t_{2}<t_{1}$$ is symmetrical w.r.t. indices $$\left(\cdot\right)_{1},\left(\cdot\right)_{2}$$. Also, the prerequisite for combining two summarizers is $$\alpha_{1}=\alpha_{2}=\alpha$$ (history decay is the same).

But enough of midnight oil burning.

Online mean summarizer for binomial distribution with irregular sampling and arbitrary bias

Suppose we have a task to provide biased estimate on probability sampled by $$\left(x_{1},...x_{n}\right),\, x_{i}\in\left\{ 0,1\right\}$$. Traditional solution for a biased estimate is mean of a conjugate prior. In this case it would be mean of the beta-distribution \[P_{n}=\mathbb{E}\left(\mu_{B}\right)=\frac{\alpha}{\alpha+\beta},\] where $$\alpha$$ equals to number of positive outcomes +1 and $$\beta$$ is taken as number of negative outcomes +1.

The obvious intention of biased estimate is to converge on certain bias $$P_{0}$$ in absence of a good sample data: \[P_{0}=\lim_{n\rightarrow0}P_{n}.\]

In case of standard beta-distribution $$P_{0}={1\over{2}}$$, which of course is not terribly useful in practice. In practice we may have much better 'guesses' for our initial estimate. So in order to allow arbitrary parameterization for $$P_{0}$$, let's denote $$\alpha=n_{+}+b_{+}$$ and $$\beta=n_{-}+b_{-}$$, $$n_{-}$$ being the number of negative observations and $$n_{+}$$ being the number of positive observations. Values $$b_{+}$$ and $$b_{-}$$ thus express our initial bias towards positive and negative outcome of final estimate, respectively.

Then our online summarizer could be presented as \[P_{n}=\frac{b_{+}+\sum_{i}^{n}x_{i}}{n+b_{-}+b_{+}}\] and it's not hard to see that we can come up with heuristics allowing arbitrary bias \[P_{0}=\lim_{n\rightarrow0}P_{n}=\frac{b_{+}}{b_{+}+b_{-}}\] while keeping $$b_{+}+b_{-}=2$$ per standard beta-distribution.

That's the biased summarizer I have been using until I saw this post.

If our sampling $$\left(x_{1},...x_{n}\right),\, x_{i}\in\left\{ 0,1\right\}$$ is also taken at times $$t_{1},...t_{n}$$ then we can factor in exponential phase-out for more distant history while taking more recent history into account. Also, instead of converging onto our bias $$P_{0}$$ when we just have a lack of history, we can also converge on it if we have lack of recent history: \[\lim_{n\rightarrow0}P=\frac{b_{+}}{b_{+}+b_{-}}=P_{0}\] and \[\lim_{\left(t_{n}-t_{n-1}\right)\rightarrow+\infty}P=\frac{b_{+}}{b_{+}+b_{-}}=P_{0}.\]

Cool. How do we exactly do that in a convenient way?

We can modify result from Ted Dunning's post I referenced above in the following way:\[P=\frac{b_{+}+\sum_{i=1}^{n}x_{i}e^{-\left(t-t_{i}\right)/\alpha}}{b_{-}+b_{+}+\sum_{i=1}^{n}e^{-\left(t-t_{i}\right)/\alpha}}.\]

It's not hard to see that our goals described by limits above would hold with this solution.

The iterative solution for that would be \[\pi_{n+1}=e^{-\left(t_{n+1}-t_{n}\right)/\alpha},\] \[w_{n+1}=1+\pi_{n+1}w_{n},\] \[s_{n+1}=x_{n+1}+\pi_{n+1}s_{n},\] \[P_{n+1}=\frac{b_{+}+s_{n+1}}{b_{+}+b_{-}+w_{n+1}}.\]

We also have to go by selecting our $$b_{+}$$ and $$b_{-}$$ values more carefully, since value of samples is exponentially decreasing. So it stands to reason we want to decrease effect of our bias based on the amount of history, exponentially weighted, as well. Suppose we have a metric $$\varepsilon$$ that denotes amount of non-significant history for the purposes of biasing. Then we want to modify condition $$b_{+}+b_{-}=2$$ based on non-weighted history by using exponential function average \[b_{+}+b_{-}=2\cdot\frac{\int_{0}^{-\ln\varepsilon}e^{-t}dt}{-\ln\varepsilon}=2\cdot\frac{\varepsilon-1}{\ln\varepsilon}.\] This gives us solution for bias parameters as \[\begin{cases}b_{+}=2P_{0}\frac{\varepsilon-1}{\ln\varepsilon};\\b_{-}=2\left(1-P_{0}\right)\frac{\varepsilon-1}{\ln\varepsilon}.\end{cases}\]

Also, having to specify $$\alpha$$ is weird. Most people like to look at it as span of useful history $$\Delta{t}$$ and margin $$m$$, in exponential scale, when the rest of history is not deemed very useful. Good default value for $$m$$ is perhaps 0.01 (1%). Then we can compute $$\alpha$$ as \[\alpha=\frac{-\Delta t}{\ln m}.\]

So, building the summarizer, we have iterative state represented by $$\left\{ w_{n},s_{n},t_{n}\right\}$$ and non-default constructor that accepts $$\Delta{t},m,P_{0}$$ and $$\varepsilon$$ and transforms them into parameters of the summarizer : $$b_{-},b_{+},\alpha$$ per above.

The summarizer needs to implement 2 methods: $$\mathrm{pnow}(t)$$ and $$\mathrm{update}(t_{n},x_{n})$$. The implementation of the update() method pretty much follows from all the above. The implementation of pnow() just needs to compute probability for current moment (we don't have an observation and may be more biased if we did not see recent observations). Method pnow() just needs to do the same estimate as update() assuming $$t_{n}$$ as 'now' and $$s_{n+1}=\pi_{n+1}s_{n}$$ and $$w_{n+1}=\pi_{n+1}w_{n}$$ and  without actually updating the state $$\left\{ w_{n},s_{n},t_{n}\right\}$$.

Did I say that some version of constructor could assume default values for $$\varepsilon$$ and $$m$$? I use $$\varepsilon=0.5$$ and $$m=0.01$$ as default. So one of good versions of constructor is the one that accepts bias aka 'null hypothesis' $$P_{0}$$, and the span of useful history $$\Delta t$$ to phase out on.

So a bit of coding, and we start converging on a pretty picture in no time.

(And yes, I do take all my pictures myself as well as doing the coding part.)

Git, Github and committing to ASF svn

There has been a discussion around ASF what the best practices of git/github/ASF svn might be. I am ready to offer my own variation.

It's been some time since ASF (Apache Software Foundation) enabled git mirrors of its SVN repository. Github has been mirroring those for some time.

E.g. Mahout project has apache git url git://git.apache.org/mahout.git which in turn is mirrored in Github as git@github.apache.org:apache/mahout.git. Being a Mahout committer myself, I will use it further as an example.

Consequently, it is possible to use Github's collaboration framework to collaborate on individual project issues, accept additions from contributors, merging/branching even more etc. In other words, enable individual jira issue to have its own commit history without really exposing this history to the final project. Another benefit is to have all the power of 3- or whatever-way rebases and merges git provides instead of having a single patch which often gets out of sync with the trunk.

And of course, ability of git to be a 'distributed' system. Your own copy is always your own play yard as well with full power of ... and so on and so forth.

So.. how do we sync apache svn with git branch?

We use git's paradigm of upstream branch, of course, with svn trunk being an upstream. Sort of.

Of course, svn doesn't support upstream tracking, not directly anyway.

First, let's clone the project and set up svn 'upstream' :

git clone git://git.apache.org/mahout.git mahout-svn
git svn init --prefix=origin/ --tags=tags --trunk=trunk --branches=branches https://svn.apache.org/repos/asf/mahout
git svn rebase

I used to run 'git svn clone...' here, but as this wonderful document points out, you can save a lot of time by cloning apache git instead of svn.

Also don't forget to install the authors file from http://git.apache.org/authors.txt by running

git config svn.authorsfile ".git/authors.txt"

assuming that's the path where you put it. Otherwise, your local git commit metadata will be different from that of created in git.apache.org and as a result, that commit's hash will also be different. "svn rebase" would reconcile it though, it seems.

We could also use 'fork' capability in github to create our repository clone, i think. It would be yet even more faster. Now that i have already a cloned repository, I cannot clone yet one more since Github doesn't seem to allow to have more than 1 fork of another repository. Oh well...

So, we are now in our 'holy cow', trunk branch. We will not work here but only use it for commits. It is also recommended to configure git-svn with svn 'authors' file per that info in apache wiki above, for committers.

Now, say we want to create a branch in Github repository, git@github.com:dlyubimov/mahout-commits which we pre-created using Github tools.

git remote add github git@github.com:dlyubimov/mahout-commits
git checkout trunk
git checkout -b MAHOUT-622  
git push -u github 

1-- creates remote corresponding to github repository, in local repository.
2-- make sure we are on local branch 'trunk'
3-- creates local branch MAHOUT-622 based on current (trunk) branch and switches to it
4-- pushes (in this case, also creates) current local MAHOUT-622 branch to github/MAHOUT-622 and also sets up upstream tracking for it. This would take significantly less time if we used Github's 'fork' as mentioned above.

Now we can create some commit history for github/MAHOUT-622 using Github collaboration (pull requests, etc.) So we are reasonably satisfied with the branch and want to commit it to trunk. Here is what we do :

git checkout MAHOUT-622
git pull
git checkout trunk
git svn rebase
git merge MAHOUT-622 --squash
git commit -m "MAHOUT-622: cleaning up dependencies"
git log -p -1 
git svn dcommit

1-- switch to MAHOUT-622
2-- fast-forward to latest changes in github remote branch
3-- switch to trunk
4-- pull latest changes from svn (just in case)
5-- merge all changes from MAHOUT-622 onto trunk tree without committing. At this point, if any merge conflicts are reported, resolve them (perhaps using 'git mergetool') if necessary;
6-- format future single svn commit for MAHOUT-622 issue as a single commit
7-- optionally check the patch of what we just changed before pushing it to svn. Also, since we just merged it here, i.e. potentially added some more changes to the original patch, do all the due diligence here: check for build, tests passing, etc. whatever else committers do.
8-- push commit to svn. European users are recommended to use --no-rebase option while doing this and then explicitly rebase several seconds later.

Alternatively, in step 5 we can merge directly from remote github branch as:

git checkout trunk
git svn rebase
git fetch github
git merge github/MAHOUT-622 --squash
git commit -m "MAHOUT-622: cleaning up dependencies"
git log -p -1 
git svn dcommit

Note that in this case, instead of pulling, we need to run 'git fetch github' in order to update remote commit tree in the local reference repo.

If we wanted to publish a patch into JIRA so that others could review it, we could do it by running 'git diff trunk MAHOUT-622' or just 'git diff trunk' if we are on MAHOUT-622 already.

Yet another use case is when we want to merge upstream svn changes into the issue we are working on :

git checkout trunk
git svn rebase
git checkout MAHOUT-622
git pull 
git merge trunk
git push

1-- switch to the 'holy cow'
2-- pull latest from 'upstream' ASF svn
3-- switch to MAHOUT-622
4-- pull latest remote collaborations to local MAHOUT-622
5-- merge trunk onto issue branch. If conflicts are reported, resolve them (perhaps using 'git mergetool');
6-- push updates to the remote github issue branch. Since we set up the issue branch to track github's branch, we can just use the default form of 'git push'.

Streaming QR decomposition for MapReduce part III: Collecting Q by Induction

In Part II, as a part of bottom-up divide-and-conquer parallelization strategy, one problem emerged that we need to solve at every node of bottom-up tree:

$$\left(\begin{matrix}\left(\begin{matrix}\mathbf{Q}_{1}\\\mathbf{Q}_{2}\\\vdots\\\mathbf{Q}\end{matrix}\right)\mathbf{R}_{1}\\\left(\begin{matrix}\mathbf{Q}\\\mathbf{Q}\\\vdots\\\mathbf{Q}\end{matrix}\right)\mathbf{R}_{2}\\\vdots\\\left(\begin{matrix}\mathbf{Q}\\\mathbf{Q}\\\vdots\\\mathbf{Q}_{z}\end{matrix}\right)\mathbf{R}_{n}\end{matrix}\right)=\left(\begin{matrix}\mathbf{\hat{Q}}\\\mathbf{\hat{Q}}\\\vdots\\\mathbf{\hat{Q}}\end{matrix}\right)\mathbf{\hat{R}}$$        (1)       

Let's simplify the problem for a moment (as it turns out, vertical blocking of Q is trivial as they are always produced by application of column-wise Givens operations which we can 'replay' on any Q block).

$$\left(\begin{matrix}\mathbf{Q}_{1}\mathbf{R}_{1}\\\mathbf{Q}_{2}\mathbf{R}_{2}\\\cdots\\\mathbf{Q}_{z}\mathbf{R}_{z}\end{matrix}\right)\Rightarrow\left(\begin{matrix}\mathbf{\hat{Q}}_{1}\\\mathbf{\hat{Q}}_{2}\\\cdots\\\mathbf{\hat{Q}}_{z}\end{matrix}\right)\hat{\mathbf{R}}$$     (2)

Also let's decompose standard Givens QR operation in to two:

1) producing Givens transformation sequence on some input $m\times n$, $m>n$  A as $$\left(\mathbf{G}_{1},\mathbf{G}_{2}...\mathbf{G}_{k}\right)=\mathrm{givens\_qr}\left(\mathbf{A}\right)$$.  I will denote product of all Givens operations obtained this way as $$\prod_{i}\mathbf{G}_{i}\equiv\prod\mathrm{givens\_qr}\left(\mathbf{A}\right)$$.

2) Applying product of Givens operations on input produces already familiar result
\[\left(\prod_{i}\mathbf{G}_{i}\right)^{\top}\mathbf{A}=\left(\begin{matrix}\mathbf{R}\\\mathbf{Z}\end{matrix}\right)\]

from which it follows that thin QR's R can be written as
\[\mathbf{R}=\left[\left(\prod_{i}\mathbf{G}_{i}\right)^{\top}\mathbf{A}\right]\left(1:n,:\right)=\left[\left(\prod\mathrm{givens\_qr}\left(\mathbf{A}\right)\right)^{\top}\mathbf{A}\right]\left(1:n,:\right)\]
using Golub/Van Loan's block notation.

(Golub/Van Loan's block notation in form of A(a₁:a₂, b₁:b₂) means "crop of a matrix A, rows a₁ to a₂ and columns b₁ to b₂". For half- and full-open intervals the bounds are just omitted).

'Thick' Q is produced by $$\mathbf{I}\left(\prod_{i}\mathbf{G}_{i}\right)$$ and hence 'thin Q' is
\[\mathbf{Q}=\left[\mathbf{I}\left(\prod_{i}\mathbf{G}_{i}\right)\right]\left(:,1:n\right)=\left[\mathbf{I}\left(\prod_{i}\mathrm{givens\_qr}\left(\mathbf{A}\right)\right)\right]\left(:,1:n\right),\mathbf{\,\, I}\in\mathbb{R}^{m\times m}.\]

Now that we laid out all notations, let's get to the gist. As I mentioned before, I continue building algorithm by induction. Let's consider case (2) for $$z=2$$. Then it can be demonstrated that the following is the solution as it is equivalent to full Givens QR with a rearranged but legitimate order of Givens operations:

$$\mathbf{\hat{R}}=\left\{\left[\prod\mathrm{givens\_qr}\left(\begin{matrix}\mathbf{R}_{1}\\\mathbf{R}_{2}\end{matrix}\right)\right]^{\top}\cdot\left(\begin{matrix}\mathbf{R}_{1}\\\mathbf{R}_{2}\end{matrix}\right)\right\} \left(1:n,:\right)$$,

$$\mathbf{\hat{Q}=}\left[\left(\begin{matrix}\mathbf{Q}_{1} & \mathbf{Z}\\\mathbf{Z} & \mathbf{Q}_{2}\end{matrix}\right)\cdot\prod_{i}\mathrm{givens\_qr}\left(\begin{matrix}\mathbf{R}_{1}\\\mathbf{R}_{2}\end{matrix}\right)\right]\left(:,1:n\right)$$.

Let's denote function of 2-block computation of $$\mathbf{\hat{R}}$$ as 

$$\mathrm{rhat}\left(\mathbf{R}_{i},\mathbf{R}_{j}\right)=\left\{\left[\prod\mathrm{givens\_qr}\left(\begin{matrix}\mathbf{R}_{i}\\\mathbf{R}_{j}\end{matrix}\right)\right]^{\top}\cdot\left(\begin{matrix}\mathbf{R}_{i}\\\mathbf{R}_{j}\end{matrix}\right)\right\} \left(1:n,:\right)$$.

Still working on the trivial case of induction: Now note that since Givens operations for calculation of $$\mathbf{\hat{Q}}$$ are applied pairwise to the columns of the accumulator matrix, i.e. independently for each row, we can split computation of $$\mathbf{\hat{Q}}$$ and apply it to any combination of vertical blocks of $$\left(\begin{matrix}\mathbf{Q}_{1} & \mathbf{Z}\\\mathbf{Z} & \mathbf{Q}_{2}\end{matrix}\right)$$. Let's say we decide to split it into 2 blocks with as many rows as in Q₁ and Q₂ and get back to block-wise formulas sought for solution of (2):

$$\mathbf{\hat{Q}_{1}=}\left[\left(\begin{matrix}\mathbf{Q}_{1} & \mathbf{Z}\end{matrix}\right)\cdot\prod_{i}\mathrm{givens\_qr}\left(\begin{matrix}\mathbf{R}_{1}\\\mathbf{R}_{2}\end{matrix}\right)\right]\left(:,1:n\right),$$

$$\mathbf{\hat{Q}_{2}=}\left[\left(\begin{matrix}\mathbf{Z} & \mathbf{Q}_{2}\end{matrix}\right)\cdot\prod_{i}\mathrm{givens\_qr}\left(\begin{matrix}\mathbf{R}_{1}\\\mathbf{R}_{2}\end{matrix}\right)\right]\left(:,1:n\right).$$      (3)

Let's denote those transformations in more general form as 

$$\mathrm{qhat\_down}\left(\mathbf{Q}_{i},\mathbf{R}_{i},\mathbf{R}_{j}\right)=\left[\left(\begin{matrix}\mathbf{Q}_{i} & \mathbf{Z}\end{matrix}\right)\cdot\prod\mathrm{givens\_qr}\left(\begin{matrix}\mathbf{R}_{i}\\\mathbf{R}_{j}\end{matrix}\right)\right]\left(:,1:n\right)$$

and

$$\mathrm{qhat\_up}\left(\mathbf{Q}_{i},\mathbf{R}_{i},\mathbf{R}_{j}\right)=\left[\left(\begin{matrix}\mathbf{Z} & \mathbf{Q}_{i}\end{matrix}\right)\cdot\prod\mathrm{givens\_qr}\left(\begin{matrix}\mathbf{R}_{j}\\\mathbf{R}_{i}\end{matrix}\right)\right]\left(:,1:n\right)$$

then we can rewrite (3) as

$$\mathbf{\hat{Q}}=\left(\begin{matrix}\mathbf{\hat{Q}}_{1}\\\mathbf{\hat{Q}}_{2}\end{matrix}\right)=\left(\begin{matrix}\mathrm{qhat\_down}\left(\mathbf{Q}_{1},\mathbf{R}_{1},\mathbf{R}_{2}\right)\\\mathrm{qhat\_up}\left(\mathbf{Q}_{2},\mathbf{R}_{2},\mathbf{R}_{1}\right)\end{matrix}\right)$$,

$$\mathbf{\hat{R}}=\mathrm{rhat}\left(\begin{matrix}\mathbf{R}_{1}\\\mathbf{R}_{2}\end{matrix}\right)$$

This is our solution for (2) of trivial case (z=2).

I will show solution for z=3 and then just will give the final solution without a proof.

Case z=3:

Note that functions qhat_up() and qhat_down() also implicitly produce intermediate rhat() products during their computation, so we want to 'enhance' them to capture that rhat() result as well. We will denote 'evolving' intermediate results Q and R as a sequence $$\left(\mathbf{\tilde{Q}},\mathbf{\tilde{R}}\right)$$:

$$\mathrm{qrhat\_down}\left(\mathbf{Q}_{i},\mathbf{R}_{i},\mathbf{R}_{j}\right)=\left(\mathbf{\tilde{Q}},\mathbf{\tilde{R}}\right)=\left(\mathrm{qhat\_down}\left(\mathbf{Q}_{i},\mathbf{R}_{i},\mathbf{R}_{j}\right),\mathrm{rhat}\left(\mathbf{R}_{i},\mathbf{R}_{j}\right)\right)$$,

$$\mathrm{qrhat\_up}\left(\mathbf{Q}_{i},\mathbf{R}_{i},\mathbf{R}_{j}\right)=\left(\mathbf{\tilde{Q}},\mathbf{\tilde{R}}\right)=\left(\mathrm{qhat\_up}\left(\mathbf{Q}_{i},\mathbf{R}_{i},\mathbf{R}_{j}\right),\mathrm{rhat}\left(\mathbf{R}_{j},\mathbf{R}_{i}\right)\right)$$.

Then, using those notations, we can write solution for (2) z=3 as

$$\mathbf{\hat{Q}}=\left(\begin{matrix}\mathbf{\hat{Q}}_{1}\\\mathbf{\hat{Q}}_{2}\\\mathbf{\hat{Q}}_{3}\end{matrix}\right)=\left(\begin{matrix}\mathrm{qrhat\_down}\left(\mathrm{qrhat\_down}\left(\mathbf{Q}_{1},\mathbf{R}_{1},\mathbf{R}_{2}\right),\mathbf{R}_{3}\right).\mathbf{\tilde{Q}}\\\mathrm{qrhat\_down}\left(\mathrm{qrhat\_up}\left(\mathbf{Q}_{2},\mathbf{R}_{2},\mathbf{R}_{1}\right),\mathbf{R}_{3}\right).\mathbf{\tilde{Q}}\\\mathrm{qrhat\_up}\left(\mathbf{Q}_{3},\mathbf{R}_{3},\mathrm{rhat}\left(\mathbf{R}_{1},\mathbf{R}_{2}\right)\right).\mathbf{\tilde{Q}}\end{matrix}\right)$$.

Note that algorithm for computing $$\mathbf{\hat{Q}}_i$$ requires input of Q_i, iterator[(R₁,R₂...R_z)].

General solution for (2) for any z is built by induction as the following algorithm:

This algorithm is still sub-efficient as for the entire matrix (2) it computes more rhat() operations than needed, and can be optimized to reduce those operations when considering final solution for (1),  but that's probably too many details for this post for now. Actual details are found in code and my working notes.

Finally, a word about how to transform solution for (2) into a solution for (1).

We can rewrite indexes in (1) as

$$\left(\begin{matrix}\left(\begin{matrix}\mathbf{Q}_{1}\\\mathbf{Q}_{2}\\\cdots\\\mathbf{Q}\end{matrix}\right)\mathbf{R}_{1}\\\left(\begin{matrix}\mathbf{Q}\\\mathbf{Q}\\\cdots\\\mathbf{Q}\end{matrix}\right)\mathbf{R}_{2}\\\cdots\\\left(\begin{matrix}\mathbf{Q}\\\mathbf{Q}\\\cdots\\\mathbf{Q}_{z}\end{matrix}\right)\mathbf{R}_{n}\end{matrix}\right)\equiv\left(\begin{matrix}\left(\begin{matrix}\mathbf{Q}_{11}\\\mathbf{Q}_{12}\\\cdots\\\mathbf{Q}_{1k}\end{matrix}\right)\mathbf{R}_{1}\\\left(\begin{matrix}\mathbf{Q}_{21}\\\mathbf{Q}_{22}\\\cdots\\\mathbf{Q}_{2k}\end{matrix}\right)\mathbf{R}_{2}\\\cdots\\\left(\begin{matrix}\mathbf{Q}_{n1}\\\mathbf{Q}_{n2}\\\cdots\\\mathbf{Q}_{nk}\end{matrix}\right)\mathbf{R}_{n}\end{matrix}\right)$$

And then regroup that into k independent tasks solving

$$\left(\begin{matrix}\mathbf{Q}_{1i}\mathbf{R}_{1}\\\mathbf{Q}_{2i}\mathbf{R}_{2}\\\cdots\\\mathbf{Q}_{zi}\mathbf{R}_{z}\end{matrix}\right)\Rightarrow\left(\begin{matrix}\mathbf{\hat{Q}}_{1i}\\\mathbf{\hat{Q}}_{2i}\\\cdots\\\mathbf{\hat{Q}}_{zi}\end{matrix}\right)\hat{\mathbf{R}}$$,

$$i\in1..k,$$

which can be solved via solution for (2).

That may also be one of  parallelization strategies.

That wasn't so hard, was it? It's possible I am performing a subpar re-tracing of some existing research or method here. But I don't think there's much work about how to fit  QR onto parallel batch machinery.

I hope that some horizons perhaps became a little clearer.

SSVD Command Line usage

Here's the doc, also attached to Mahout-593. At some point wiki update is due. When we know what it is all good for.



Also, from my email regarding -s parameter:

There are 2 cases where you might want to adjust -s :

1 -- if you running really huge input that produces more than 1000 or
so map tasks and/or that is causing OOM in some tasks in some
situations. It looks like your input is far from that now.

2 -- if you have quite wide input -- realistically more than 30k
non-zero elements in a row. The way current algorithm works, it tries
to do blocking QR of stochastically projected rows in the mappers
which means it needs to read at least k+p rows in each split (map
task). This can be fixed and i have a branch that should eventually
address this.  In your case, if there happen to be splits that contain
less than 110 rows of input, that would be the case where you might
want to start setting -s greater than DFS block size (64mb) but it has
no effect if it's less than that (which is why hadoop calls it
_minimum_ split size). I don't remember hadoop's definition of this
parameter, i think it is in bytes, so that means you probably need to
specify something like 100,000,000 to start seeing decrease in number
of the map tasks. But honestly i never tried this yet since i never
had input wide enough to require this.

Streaming QR decomposition for MapReduce part II: Bottom-up divide-and-conquer overview

In Part I of "Streaming QR decomposition for MapReduce" I touched a little bit about how bottom-up collection of Q is done. I think I would like to add a couple of figures in attempt to give a little more details and clarifications how it is done.

Collecting Q: Divide-And-Conquer: a tad more details


As I mentioned before, outer QR step is essentially bottom-up n-indegree Divide-And-Conquer algorithm.

\[\mathbf{Y}=\left(\begin{matrix}\cdots\\\mathbf{Y}_{i}\\\mathbf{Y}_{i+1}\\\mathbf{Y}_{i+2}\\\cdots\end{matrix}\right)=\begin{matrix}\cdots\\\left.\begin{matrix}\mathbf{Q} & \mathbf{R}\\\mathbf{Q} & \mathbf{R}\\\cdots & \cdots\\\mathbf{Q} & \mathbf{R}\end{matrix}\right]\\\\\left.\begin{matrix}\mathbf{Q} & \mathbf{R}\\\mathbf{Q} & \mathbf{R}\\\cdots & \cdots\\\mathbf{Q} & \mathbf{R}\end{matrix}\right]\\\\\left.\begin{matrix}\mathbf{Q} & \mathbf{R}\\\mathbf{Q} & \mathbf{R}\\\cdots & \cdots\\\mathbf{Q} & \mathbf{R}\end{matrix}\right]\\\cdots\end{matrix}\Rightarrow\begin{matrix}\cdots\\\left.\begin{matrix}\cdots & \mathbf{\cdots}\\\left.\begin{matrix}\mathbf{Q}\\\mathbf{Q}\\\cdots\\\mathbf{Q}\end{matrix}\right] & \mathbf{R}\\\\\left.\begin{matrix}\mathbf{Q}\\\mathbf{Q}\\\cdots\\\mathbf{Q}\end{matrix}\right] & \mathbf{R}\\\\\left.\begin{matrix}\mathbf{Q}\\\mathbf{Q}\\\cdots\\\mathbf{Q}\end{matrix}\right] & \mathbf{R}\\\cdots & \cdots\end{matrix}\right]\\\cdots\end{matrix}\Rightarrow\cdots\Rightarrow\begin{matrix}\left.\begin{matrix}\mathbf{\cdots}\\\mathbf{\hat{Q}}\\\mathbf{\hat{Q}}\\\cdots\\\mathbf{\hat{Q}}\\\\\mathbf{\hat{Q}}\\\mathbf{\hat{Q}}\\\cdots\\\mathbf{\hat{Q}}\\\\\mathbf{\hat{Q}}\\\mathbf{\hat{Q}}\\\cdots\\\mathbf{\hat{Q}}\\\cdots\end{matrix}\right] & \hat{\mathbf{R}}\end{matrix}\]

 The way generalized Givens thin QR works, it applies a number of Givens transformations on a Q-accumulator matrix given some initial tall matrix (in this case, Y blocks) until the tall matrix is reduced to a form of $\left(\begin{matrix}\mathbf{R}\\\mathbf{Z}\end{matrix}\right)$ , where Z has all zero elements and R is an upper-triangular.   Initial value for Q-accumulator matrix is chosen as I (square identity matrix). To tranform result to a 'thin' QR result, only first n columns of accumulator are taken (that becomes Q of the decomposition) and R part is taken as the second matrix of the thin decomposition.

In our case, in order to enable Givens recursive use, we generalize standard Givens algorithm by making it  accept two parameters: pre-existing 'thick qr' Q accumulator instead of I (in reality only shifting 'thin' Q accumulator is sufficient) and the tall input matrix as second input.

Bottom-up divide-and-conquer algorithm hence can be described as follows:

• First, we split all input (Y) into vertical blocks (see illustration above). For each block, we compute slightly modified version of standard Givens QR which is optimized for sequential (streaming) access to the tall input block elements without running into big memory requirements. The output of this step is bunch of (Q,R) pairs (each pair corresponds to the original block). For benefit of subsequent step, we will consider this block data as ({Q}, R) where sequence of Q blocks is denoted as {Q} and contains only one matrix. (Actually math stuff normally uses parenthesis () to distinguish sequences from sets, but I feel using {} notation is much more expressive in this case).
  
  
• Second, we group the input of form {Q}, R into new groups, each group thus would be denoted as {{Q},R} such that number of R in each group doesn't exceed certain maximum bottom-up indegree limit N (usually 1000 for 1G RAM solvers).
  
  
• Third, for each group formed in step 2 above, we run row-wise Givens solver, which produces ("merges")  a valid Givens solution for  each group {{Q},R} → {Q},R. Essentually this solves "QR reduction" at nodes of divide-and-conquer:
  
  
  ,
  which is essentially equivalent to solving a number of individual block-wise QR decompositions over a one ore more initial blocks of Y into single one
  .
  This algorithm is called 'compute QHatSequence2' (I think) in my detailed notes and is the one that builds by induction. Hat sign is to denote final QR blocks. I plan to discuss details of that algorithm in "part III" of this blog.
  
  
• Fourth, we consider each result of step 3 to form a sequence again and restart form step 2 and repeat it from there unless the number of groups is 1, which would produce our final Q blocks and single R as a result.

Parallelization Strategy

1. We can consider the iterations above from the point of view of individual block Y_i as a set of isolated parallel steps. We can view the entire computation as a series of independent computations over Q_i,{R₁...R_n}. It also turns out we always consume the sequence of R in the same order and never have to load more than one upper triangular matrix in the memory, so the parameters of such algorithms can actually be [Q_i, iterator{R₁...R_n}]. The algorithm produces new version of Q block, and the last solver in a group would produce a new R item for the next step R-sequence (as can be demonstrated further on). Then solver is reloaded with next R sequence and runs again (until we are left with just one R). That would be next map-only job run but each run reduces total number of Rs thousands of times. so 3 map-only runs can handle 1 billion blocks (besides the last run is really combined with the next step, computation of $\mathbf{B}^{\top}$, so we save at least one iteration setup overhead here.
2. Another small enhancement is that it turns out each Q computation doesn't need entire {R₁...R_n}. Let's rewrite Q block indexes into 2-variable index to reflect their group and subgroup indices:\[\left(\begin{matrix}\left(\begin{matrix}\mathbf{Q}_{1}\\\mathbf{Q}_{2}\\\cdots\\\mathbf{Q}\end{matrix}\right)\mathbf{R}_{1}\\\left(\begin{matrix}\mathbf{Q}\\\mathbf{Q}\\\cdots\\\mathbf{Q}\end{matrix}\right)\mathbf{R}_{2}\\\cdots\\\left(\begin{matrix}\mathbf{Q}\\\mathbf{Q}\\\cdots\\\mathbf{Q}_{z}\end{matrix}\right)\mathbf{R}_{n}\end{matrix}\right)\equiv\left(\begin{matrix}\left(\begin{matrix}\mathbf{Q}_{11}\\\mathbf{Q}_{12}\\\cdots\\\mathbf{Q}_{1k}\end{matrix}\right)\mathbf{R}_{1}\\\left(\begin{matrix}\mathbf{Q}_{21}\\\mathbf{Q}_{22}\\\cdots\\\mathbf{Q}_{2k}\end{matrix}\right)\mathbf{R}_{2}\\\cdots\\\left(\begin{matrix}\mathbf{Q}_{n1}\\\mathbf{Q}_{n2}\\\cdots\\\mathbf{Q}_{nk}\end{matrix}\right)\mathbf{R}_{n}\end{matrix}\right).\] It turns out that computation over Q_1i and Q_2i  requires entire {R₁...R_n} but computation over Q_3i requires $\left\{ \mathrm{GivensQR}\left[\left(\begin{matrix}\mathbf{R}_{1}\\\mathbf{R}_{2}\end{matrix}\right)\right].\mathbf{R},\mathbf{R}_{3}...\mathbf{R}_{n}\right\}$, i.e only n-1 upper-triangulars. Going on, сomputation over Q_4i requires sequence $\left\{ \mathrm{GivensQR}\left[\left(\begin{matrix}\mathrm{GivensQR}\left[\left(\begin{matrix}\mathbf{R}_{1}\\\mathbf{R}_{2}\end{matrix}\right)\right].\mathbf{R}\\\mathbf{R}_{3}\end{matrix}\right)\right].\mathbf{R},\mathbf{R}_{4}...\mathbf{R}_{n}\right\}$. And so on, with the Q_n requiring only two upper-triangular arguments. That means that we quite legally can split computation into at most k independent jobs, each i-th parallel job computing blocks Q_1i, Q_2i...Q_ni in sequence while also reducing R sequence per above. To aid parallelization even more, we can actually choose k to be whatever we want: inside the algorithm, matrix $\left(\begin{matrix}\mathbf{Q}_{i1}\\\mathbf{Q}_{i2}\\\cdots\\\mathbf{Q}_{ik}\end{matrix}\right)$ is only ever transformed by applying column-wise Givens operations and horizontal shifts. Hence, it doesn't matter how it is being split into blocks, we can regroup rows there in any way to come up with a degree of parallelism k we are comfortable with.

Looks pretty complicated, huh? I guess it's not most complicated part yet though. Besides, I haven't implemented bottom-up approach 100% yet. I only implemented 2-step hierarchy (i.e. we can have n×n Y blocks initially only). I think it is kind of enough for what my company does, but it can be worked on to insert more steps in between to enable as many blocks as we want.

We've just started on the path uphill.