Full code:
http://theparticle.com/applets/ml/least_squares_nonlinear/pMatrix.java

    /**
     * non-negative matrix factorization
     * 
     * w.mult(h) = this
     * w is weights matrix.
     * h is feature matrix.
     *
     * k = size of factor matrices.
     * n is NUMBER OF ITERATIONS.
     * e is ERROR (once we're below that, we can return)
     * 
     * returns w,h.
     */
    public pMatrix[] nmf(int k,int n,double e){
        pMatrix w = urandom(rows,k,1,n);
        pMatrix h = urandom(k,cols,1,n);
        for(int i=0;i<n;i++){
            // compute output.
            pMatrix wh = w.mult(h);
            double cost = distSq(wh);

            // if found solution
            if(cost < e)
                break;

            // update feature matrix.
            pMatrix wt = w.T();
            pMatrix hn = wt.mult(this);
            pMatrix hd = wt.mult(wh);
            h.smultEq(hn.smultEq(hd.sinvEq()));

            // update weights matrix
            pMatrix ht = h.T();
            pMatrix wn = mult(ht);
            pMatrix wd = w.mult(h).mult(ht);
            w.smultEq(wn.smultEq(wd.sinvEq()));
        }
        return new pMatrix[]{ w, h };
    }