libSkylark documentation

Quick Start Guide

«  Introduction   ::   Contents   ::   Building from Source  »

Quick Start Guide

There are two options for quickly installing libskylark in Linux-x86_64 systems (i.e. Linux distribution independent):

Installing libskylark from conda channels

Open a terminal and type (copy-paste):

# Get miniconda for python 2.7
curl -L -O https://repo.continuum.io/miniconda/Miniconda2-latest-Linux-x86_64.sh
bash Miniconda2-latest-Linux-x86_64.sh -b -p ${HOME}/miniconda2

# Get the libskylark channel
curl -L -O https://github.com/xdata-skylark/libskylark/releases/download/v0.20/channel.tar.gz
cd ${HOME}; tar -xvzf channel.tar.gz

# Install fftw and libskylark conda packages
export PATH=${HOME}/miniconda2/bin:${PATH}
conda install fftw --yes --channel conda-forge
conda install libskylark --yes --channel file://${HOME}/channel

Trying the installation

Start an ipython shell:

ipython

and then input the following python snippet:

import El
from skylark import sketch
A = El.DistMatrix()
El.Uniform(A, 10, 10)
S = sketch.FJLT(10, 4)
SA = S * A
print [[SA.Get(i, j) for j in range(10)] for i in range(4)]

Congratulations! You have just sketched a 10x10 distributed matrix of uniform distribution random entries into a 4x10 sketched matrix by making use of the Fast Johnson-Lindenstrauss Transform (FJLT).

Now consider hard-wiring the installation path by appending this line

export PATH=${HOME}/miniconda2/bin:${PATH}

to your ~/.bashrc.

Note

Consider expanding LD_LIBRARY_PATH if intending to use example libskylark executables: export LD_LIBRARY_PATH=${HOME}/miniconda2/lib64:${HOME}/miniconda2/lib:${LD_LIBRARY_PATH}"

Installing libskylark from a single-file installer

Open a terminal and type (copy-paste):

# Download the installer and install libskylark
curl -L -O https://github.com/xdata-skylark/libskylark/releases/download/v0.20/install.sh
bash install.sh -b -p ${HOME}/libskylark
export PATH=${HOME}/libskylark/bin:${PATH}

Trying the installation

Start an ipython shell:

ipython

and then input the following python snippet:

import El
from skylark import sketch
A = El.Matrix()
El.Uniform(A, 10, 10)
S = sketch.JLT(10, 4)
SA = S * A
print SA.ToNumPy()

Congratulations! You have just sketched a 10x10 matrix of uniform distribution random entries into a 4x10 sketched matrix by making use of the Johnson-Lindenstrauss Transform (JLT).

Now consider hard-wiring the installation path by appending this line

export PATH=${HOME}/libskylark/bin:${PATH}

to your ~/.bashrc.

Note

Consider expanding LD_LIBRARY_PATH if intending to use example libskylark executables: export LD_LIBRARY_PATH=${HOME}/libskylark/lib64:${HOME}/libskylark/lib:${LD_LIBRARY_PATH}"

Vagrant

Setting up a libSkylark environment is facilitated using Vagrant and VirtualBox.

Note

Make sure to use Vagrant >= 1.6.3! The configuration file will not be compatible with older versions.

Simply follow the instructions shown below.

  1. Install VirtualBox (download VirtualBox)

  2. Install Vagrant (download Vagrant).

  3. Execute the following commands in a terminal:

    git clone https://github.com/xdata-skylark/libskylark.git
cd libskylark/vagrant/trusty64
vagrant up
vagrant ssh

Note

The vagrant up command is downloading, building and installing libsklylark and its dependencies. This step will most likely take a while to complete.

To exit the virtual environment type exit in the shell followed by vagrant halt in order to power down the virtual machine. In order to reconnect to the virtual machine, either use the VirtualBox interface to start the virtual machine and then login with the vagrant/vagrant user/password or execute

vagrant up --no-provision
vagrant ssh

Note

Depending on your Vagrant version, omitting the --no-provision argument (after a vagrant halt) is considerably slower because the vagrant up reruns parts of the bootstrap script.

Note

The following environment variables are available after login, pointing to the libSkylark source and install directory:

SKYLARK_SRC_DIR=/home/vagrant/libskylark
SKYLARK_INSTALL_DIR=/home/vagrant/install

In the next section we provide more in depth details on how the Vagrant setup works and some trouble shooting. Additionally we provide instruction on how to install on a cluster (see Cluster of vagrant-controlled VMs) and running on AWS (see Running Vagrant on AWS).

Finally, some examples are provided in the Examples of Library Usage section.

Note

You might get a warning concerning AVX instructions when running the examples

OpenBLAS : Your OS does not support AVX instructions. OpenBLAS is using Nehalem kernels as a fallback, which may give poorer performance.

because some VM providers do not support AVX instructions yet (e.g. see: https://www.virtualbox.org/ticket/11347).

Under vagrant/trusty64 we provide an easy way to materialize an 64-bit Ubuntu-based virtual machine (release 14.04 LTS, aka Trusty Tahr) with all software dependencies required for a fresh build of libSkylark.

The Vagrantfile is the input to Vagrant, a program that facilitates the management of virtual machines. Using VirtualBox as the backend virtualization software (provider) we can launch a libSkylark environment independent of the host operating system (Windows, Mac OS X and Linux host platform require exactly the same commands).

After installing VirtualBox and Vagrant the following commands suffice

git clone https://github.com/xdata-skylark/libskylark.git
cd vagrant/trusty64
vagrant up
vagrant ssh

to get a virtual machine (VM) with libSkylark and all its dependencies. If you want to customize the installed libSkylark follow the instructions (see Building libSkylark) for building and installing libSkylark.

Note

To customize the amount of memory (4096) and CPU (75%) cap that will be used for the VM, check the Vagrant file.

The vagrant box can be stopped by issuing vagrant halt.

Note

We use clean vagrant boxes, automatically downloaded from the following urls:

Note

For a Windows XP host, PuTTy ssh client should use the default connection settings: IP= 127.0.0.1, port= 2222, username= vagrant, password= vagrant.

Hint

If you get a timeout error (waiting for the virtual machine to boot)

[...]
default: Warning: Connection timeout. Retrying...
default: Warning: Connection timeout. Retrying...
Timed out while waiting for the machine to boot. This means that
Vagrant was unable to communicate with the guest machine within
the configured ("config.vm.boot_timeout" value) time period.

during vagrant up, halt the virtual machine and vagrant up again.

Hint

In case Vagrant crashes in the bootstrap phase, the safest thing is to delete the virtual machine (halt than delete the machine in the VirtualBox GUI) and start with a fresh Vagrant build.

Command-line Usage

While libSkylark is designed to be a library, it also provides access to many of its feature using standalone applications. For all applications, running with the –help options reveals the command-line options.

The following example sessions assume the following commands have been performed in advance (it downloads and extracts some demo data):

wget  https://www.dropbox.com/s/0ilbyeem2ihmvzw/data.tar.gz
tar -xvzf data.tar.gz

Approximate Singular Value Decomposition

Building libSkylark creates an executable called skylark_svd in $SKYLARK_INSTALL_DIR/bin. This executable can be used in standalone mode to compute an approximate SVD with 10 leading singular vectors as follows:

skylark_svd -k 10 --prefix usps data/usps.train

The files usps.U.txt, usps.S.txt and usps.V.txt contain the approximate SVD that was computed.

Linear Least-Squares

Building libSkylark creates an executable called skylark_linear in $SKYLARK_INSTALL_DIR/bin. This executable can be used in standalone mode to approximately solve linear least squares problems. Here is an example command-line:

skylark_linear data/cpu.train cpu.sol

The file cpu.sol.txt contains the computed approximate solution.

Note

Currently, only overdetermined least squares of dense matrices is supported. No regularization options are available. This will be relaxed in the future.

Learning Non-Linear Models

libSkylark offers two mechanisms to learn non-linear models: ADMM-based solver and kernel ridge regression solver.

Kernel Ridge Regression solver

In $SKYLARK_SRC_DIR/ml/skylark_krr.cpp, kernel ridge regression is setup to train and predict with a kernel based model for nonlinear classification and regression problems. The code incorporates various algorithms, with varying trade-offs in speed, memory and accuracy.

Building libSkylark creates an executable called skylark_krrl in $SKYLARK_INSTALL_DIR/bin. This executable can be used in standalone mode as follows.

  1. Train a classifier using the Gaussian Kernel
mpiexec -np 4 ./skylark_krr -a 1 -k 0 -g 10 -f 1000 --model model data/usps.train
  1. Test accuracy of the generated model and generate predictions
mpiexec -np 4 ./skylark_krr --outputfile output --model model usps.test

ADMM-based solver

In $SKYLARK_SRC_DIR/ml/skylark_ml.cpp, an ADMM-based solver is setup to train and predict with a randomized kernel based model for nonlinear classification and regression problems.

Building libSkylark creates an executable called skylark_ml in $SKYLARK_INSTALL_DIR/bin. This executable can be used in standalone mode as follows.

  1. Train an SVM with Randomized Gaussian Kernel
mpiexec -np 4 ./skylark_ml -g 10 -k 1 -l 2 -i 30 -f 1000 --trainfile data/usps.train --valfile data/usps.test --modelfile model
  1. Test accuracy of the generated model and generate predictions
mpiexec -np 4 ./skylark_ml --outputfile output --testfile data/usps.test --modelfile model

An output file named output.txt with the predicted labels is created.

  1. In the above, the entire test data is loaded to memory and the result is computed. This is the fast, but it is limited to cases where the test data is fixed. skylark_ml also supports an interactive (streaming) mode in which the software prompts for data line by line from stdin and outputs in prediction after each line is received. The mode is invoked by not passing either training data or testing data. An example for using this mode is the following:
cut -d ' ' -f 1 --complement data/usps.test | skylark_ml --modelfile model

Note

Interactive mode is much slower, since there is significant overhead that is repeated for each new line of data. In the future, batching will be supported, to avoid some of this overhead.

  1. It is also possible to load the model in Python and do predictions. Here is an example:
import skylark.io, skylark.ml.modeling
model = skylark.ml.modeling.LinearizedKernelModel('model')
testfile = skylark.io.libsvm('data/usps.test')
Xt, Yt = testfile.read(model.get_input_dimension())
Yp = model.predict(Xt)

Community Detection Using Seed Node

Building libSkylark creates an executable called skylark_community in $SKYLARK_INSTALL_DIR/bin. This executable can be used in standalone mode as follows. (Note that an interactive mode is also present.)

1. Prepare the input. Basically, the graph is described in a text file, with each row an edge. Each row has two strings that identify the node. The identifier can be arbitrary (does not have to be numeric). For example:

cat data/two_triangles
A1 A2
A2 A3
A1 A3
B1 B2
B2 B3
B1 B3
A1 B1
  1. Detect the community with A1 as the seed:
skylark_community --graphfile data/two_triangles --seed A1

3. The community detection algorithm is local, so it operates only on a few nodes. So, for large graphs the majority of the time for detecting a single community will be spent on just loading the file. To detect multiple communities with the graph already loaded in memory, an interactive quiet mode is provided. Here is an example:

cat data/seeds | skylark_community --graphfile data/two_triangles -i -q

Approximate Adjacency Spectral Embedding

Building libSkylark creates an executable called skylark_graph_se in $SKYLARK_INSTALL_DIR/bin. This executable can be used in standalone mode as follows.

1. Prepare the input. Basically, the graph is described in a text file, with each row an edge. Each row has two strings that identify the node. The identifier can be arbitrary (does not have to be numeric). For example:

cat data/two_triangles
A1 A2
A2 A3
A1 A3
B1 B2
B2 B3
B1 B3
A1 B1
  1. Compute a two dimensional embedding.
skylark_graph_se -k 2 data/two_triangles

3. Two files are created: out.index.txt and out.vec.txt. The vec file has the embedding vectors: each row corresponds to an index. The index file maps row index to vertices. The prefix can be called using the prefix command line option (default is out).

Cluster of vagrant-controlled VMs

Here is a simple approach for building a cluster of vagrant-controlled VMs; it works in a local setting but since it uses the bridged mode it should work in a real cluster environment as well. Please refer to this blog entry for more information on Networking in VirtualBox.

  • In Vagrantfile set:
config.vm.network :public_network

  • During vagrant up choose the ethernet interface to use.

  • After vagrant ssh do ifconfig, get the interface name - say

    eth1 - and then for the VMs at nodes 1, 2,... do:

sudo /sbin/ifconfig eth1:vm 192.168.100.1
sudo /sbin/ifconfig eth1:vm 192.168.100.2
...
  • 192.168.100.xxx will be the names to put in the “hosts” file for MPI

    daemons to use; as previously noted: vagrant/vagrant is the default user/password combination for ssh.

Running Vagrant on AWS

This follows the instructions found at:

Using the provided bootstrap.sh file in combination with the aws plugin (see link above) one can deploy on AWS. First adapt the Vagrant file as shown below:

Vagrant.configure("2") do |config|
    config.vm.box = "skylark"
    config.vm.provision :shell, :path => "bootstrap.sh"

    config.vm.provider :aws do |aws, override|
        aws.access_key_id     = "XYZ"
        aws.secret_access_key = "AAA"
        aws.keypair_name      = "keynam"
        aws.security_groups   = ["ssh-rule"]

        aws.ami           = "ami-fa9cf1ca"
        aws.region        = "us-west-2"
        aws.instance_type = "t1.micro"

        override.ssh.username         = "ubuntu"
        override.ssh.private_key_path = "aws.pem"
    end
end

Then execute:

  • vagrant plugin install vagrant-aws
  • vagrant box add skylark https://github.com/mitchellh/vagrant-aws/raw/master/dummy.box
  • vagrant up --provider=aws
  • vagrant ssh

It will take a while to compile and install everything specified in the bootstrap.sh script.

Running MPI on Hadoop/Yarn

libSkylark can be run on on a Yarn scheduled cluster using the blueprint outlined below.

Note

The installation requires admin privileges. Also make sure to satisfy the recommended dependencies.

1. Building MPICH2

Download the MPICH version specified in the dependencies (MPICH-3.1.2).

sh autogen.sh
./configure --prefix=$HOME/mpi/bin
make
make install
sudo ln -s /home/ubuntu/mpi/bin/bin/* /usr/local/bin

and copy MPI binaries to all nodes, i.e.

scp /home/ubuntu/mpi/bin/bin/* root@hadoopX:/usr/local/bin

Make sure that all the MPI executables are accessible on all the nodes.

2. Building mpich2-yarn

Finally we have all the dependencies ready and can continue to mpich2-yarn https://github.com/alibaba/mpich2-yarn:

cd ~/mpi
git clone https://github.com/alibaba/mpich2-yarn
cd mpich2-yarn
mvn clean package -Dmaven.test.skip=true

The configuration for Hadoop/Yarn can be found on on the Github project page. Note that the HDFS permission must be set correctly for the temporary folders used to distribute MPI executables.

Finally, an MPI job can be submitted, i.e. with

sudo -u yarn hadoop jar target/mpich2-yarn-1.0-SNAPSHOT.jar -a hellow -M 1024 -m 1024 -n 1

where we chose to run as the yarn user (double check permission for any user that runs mpich2-yarn). More examples can be found on the mpich2-yarn Github page.

Examples of Library Usage

Here, we provide a flavor of the library and its usage. We assume that the reader is familiar with sketching and its applications in randomized Numerical Linear algebra and Machine Learning. If not, we refer the reader to subsequent sections which provide the necessary background and references.

When libSkylark is built, executables instantiating these examples can be found under $SKYLARK_INSTALL_DIR/bin/skylark_examples and $SKYLARK_INSTALL_DIR/bin.

In addition Vagrant provisions a IPyhton notebook <http://ipython.org/notebook.html> server with access to the libSkylark Python bindings. To play with libSkylark just open your browser and direct it to localhost:6868.

Sketching

In $SKYLARK_SRC_DIR/examples/elemental.cpp an example is provided that illustrates sketching. Below, three types of sketches are illustrated in the highlighted lines:

  1. a 1D row-distributed Elemental dense matrix is sketched using the JLT (Johnson-Lindenstrauss) Transform into a local matrix
  2. a row-distributed Elemental dense matrix is sketched using FJLT (Fast Johnson-Lindenstrauss transform) involving faster FFT-like operations.
  3. a 2D block-cyclic distributed Combinatorial BLAS sparse matrix is sketched to a dense local Elemental matrix using Clarkson-Woodruff transform which involves hashing the rows.
  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
#include <iostream>
#include <functional>
#include <cstring>
#include <vector>
#include <utility>
#include <unordered_map>

#include <El.hpp>
#include <boost/mpi.hpp>

#define SKYLARK_NO_ANY
#include <skylark.hpp>

#include "utilities.hpp"
#include "parser.hpp"

/*******************************************/
namespace bmpi =  boost::mpi;
namespace skys =  skylark::sketch;
namespace skyb =  skylark::base;
/*******************************************/

/* These were declared as extern in utilities.hpp --- defining it here */
int int_params [NUM_INT_PARAMETERS];
char* chr_params[NUM_CHR_PARAMETERS];

/** Typedef DistMatrix and Matrix */
typedef std::vector<int> IntContainer;
typedef std::vector<double> DblContainer;
typedef El::DistMatrix<double, El::CIRC, El::CIRC> MatrixType;
typedef El::DistMatrix<double, El::VC, El::STAR> DistMatrixType;

int main (int argc, char** argv) {
    /* Initialize Elemental */
    El::Initialize (argc, argv);

    /* MPI sends argc and argv everywhere --- parse everywhere */
    parse_parameters (argc,argv);

    /* Initialize skylark */
    skyb::context_t context (int_params[RAND_SEED_INDEX]);

    /* Create matrices A and B */
    bmpi::communicator world;
    MPI_Comm mpi_world(world);
    El::Grid grid (mpi_world);
    El::DistMatrix<double, El::VR, El::STAR> A(grid);
    El::DistMatrix<double, El::VR, El::STAR> B(grid);

    /** Only randomization is supported for now */
    if (0==int_params[USE_RANDOM_INDEX]) {
        /** TODO: Read the entries! */
        std::cout << "We don't support reading --- yet --" << std::endl;
    } else {
        El::Uniform (A, int_params[M_INDEX], int_params[N_INDEX]);
        El::Uniform (B, int_params[M_INDEX], int_params[N_RHS_INDEX]);
    }

    /**
     * Depending on which sketch is requested, do the sketching.
     */
    if (0==strcmp("JLT", chr_params[TRANSFORM_INDEX]) ) {

        if (SKETCH_LEFT == int_params[SKETCH_DIRECTION_INDEX]) {

            /* 1. Create the sketching matrix */
            skys::JLT_t<DistMatrixType, MatrixType> JLT (int_params[M_INDEX],
                int_params[S_INDEX], context);

            /* 2. Create space for the sketched matrix */
            MatrixType sketch_A(int_params[S_INDEX], int_params[N_INDEX]);

            /* 3. Apply the transform */
            try {
                JLT.apply (A, sketch_A, skys::columnwise_tag());
            } catch (skylark::base::skylark_exception ex) {
                SKYLARK_PRINT_EXCEPTION_DETAILS(ex);
                SKYLARK_PRINT_EXCEPTION_TRACE(ex);
                errno = *(boost::get_error_info<skylark::base::error_code>(ex));
                std::cout << "Caught exception, exiting with error " << errno << std::endl;
                std::cout << skylark_strerror(errno) << std::endl;
                return errno;
            }

            /* 4. Print and see the result (if small enough) */
            if (int_params[S_INDEX] * int_params[N_INDEX] < 100 &&
                world.rank() == 0)
                El::Display(sketch_A);

            /** TODO: Do that same to B, and solve the system! */

        } else {
            std::cout << "We only have left sketching. Please retry" << std::endl;
        }
    } else if (0==strcmp("FJLT", chr_params[TRANSFORM_INDEX]) ) {
        if (SKETCH_LEFT == int_params[SKETCH_DIRECTION_INDEX]) {
            /* 1. Create the sketching matrix */
            skys::FJLT_t<DistMatrixType, MatrixType> FJLT (int_params[M_INDEX],
                int_params[S_INDEX], context);

            /* 2. Create space for the sketched matrix */
            MatrixType sketch_A(int_params[S_INDEX], int_params[N_INDEX]);

            /* 3. Apply the transform */
            FJLT.apply (A, sketch_A, skys::columnwise_tag());

            /* 4. Print and see the result */
            if (int_params[S_INDEX] * int_params[M_INDEX] < 100 &&
                world.rank() == 0)
                El::Display(sketch_A);

            /** TODO: Do that same to B, and solve the system! */


        }
    } else if (0==strcmp("CWT", chr_params[TRANSFORM_INDEX]) ) {
        if (SKETCH_LEFT == int_params[SKETCH_DIRECTION_INDEX]) {

            /* 1. Create the sketching matrix */
            skys::CWT_t<DistMatrixType, MatrixType>
                Sparse (int_params[M_INDEX], int_params[S_INDEX], context);

            /* 2. Create space for the sketched matrix */
            MatrixType sketch_A(int_params[S_INDEX], int_params[N_INDEX]);

            /* 3. Apply the transform */
            SKYLARK_BEGIN_TRY()
                Sparse.apply (A, sketch_A, skys::columnwise_tag());
            SKYLARK_END_TRY()
            SKYLARK_CATCH_AND_RETURN_ERROR_CODE();

            /* 4. Print and see the result */
            if (int_params[S_INDEX] * int_params[M_INDEX] < 100 &&
                world.rank() == 0)
                El::Display(sketch_A);

            /** TODO: Do that same to B, and solve the system! */


        }
    } else {
        std::cout << "We only have JLT/FJLT/Sparse sketching. Please retry" <<
            std::endl;
    }

    El::Finalize();

    return 0;

}

For sketching Elemental and CombBLAS matrices via Python bindings, run, for example, the following:

mpiexec -np 4 python $SKYLARK_SRC_DIR/python-skylark/skylark/examples/example_sketch.py
mpiexec -np 4 python $SKYLARK_SRC_DIR/python-skylark/skylark/examples/example_sparse_sketch.py

Note

Currently, we have stable python bindings only for sketching layer. This too might change with newer versions of Elemental.

Fast Least Square Regression

In $SKYLARK_SRC_DIR/examples/least_squares.cpp examples are provided on how least squares problems can be solved faster using sketching. One approach is of the flavor of sketch-and-solve geared towards lower-precision solutions, while the other approach uses sketching to construct a preconditioner to obtain high-precision solutions faster.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
#include <iostream>

#include <El.hpp>
#include <boost/mpi.hpp>
#include <boost/format.hpp>

#define SKYLARK_NO_ANY
#include <skylark.hpp>

const int m = 50000;
const int n = 500;


template<typename MatrixType, typename RhsType, typename SolType>
void check_solution(const MatrixType &A, const RhsType &b, const SolType &x, 
    const RhsType &r0,
    double &res, double &resAtr, double &resFac) {
    RhsType r(b);
    skylark::base::Gemv(El::NORMAL, -1.0, A, x, 1.0, r);
    res = skylark::base::Nrm2(r);

    SolType Atr(x.Height(), x.Width(), x.Grid());
    skylark::base::Gemv(El::TRANSPOSE, 1.0, A, r, 0.0, Atr);
    resAtr = skylark::base::Nrm2(Atr);

    skylark::base::Axpy(-1.0, r0, r);
    RhsType dr(b);
    skylark::base::Axpy(-1.0, r0, dr);
    resFac = skylark::base::Nrm2(r) / skylark::base::Nrm2(dr);
}

template<typename MatrixType, typename RhsType, typename SolType>
void experiment() {
    typedef MatrixType matrix_type;
    typedef RhsType rhs_type;
    typedef SolType sol_type;

    double res, resAtr, resFac;

    boost::mpi::communicator world;
    int rank = world.rank();

    skylark::base::context_t context(23234);

    // Setup problem and righthand side
    // Using Skylark's uniform generator (as opposed to Elemental's)
    // will insure the same A and b are generated regardless of the number
    // of processors.
    matrix_type A, b;
    skylark::base::UniformMatrix(A, m, n, context);
    skylark::base::UniformMatrix(b, m, 1, context);

    sol_type x(n,1);
    rhs_type r(b);

    boost::mpi::timer timer;
    double telp;

    // Solve using Elemental. Note: Elemental only supports [MC,MR]...
    El::DistMatrix<double> A1 = A, b1 = b, x1;
    timer.restart();
    El::LeastSquares(El::NORMAL, A1, b1, x1);
    telp = timer.elapsed();
    x = x1;
    check_solution(A, b, x, r, res, resAtr, resFac);
    if (rank == 0)
        std::cout << "Elemental:\t\t\t||r||_2 =  "
                  << boost::format("%.2f") % res
                  << "\t\t\t\t\t\t\t||A' * r||_2 = " << boost::format("%.2e") % resAtr
                  << "\t\tTime: " << boost::format("%.2e") % telp << " sec"
                  << std::endl;
    double res_opt = res;

    // The following computes the optimal residual (r^\star in the logs)
    skylark::base::Gemv(El::NORMAL, -1.0, A, x, 1.0, r);

#if SKYLARK_HAVE_FFTW || SKYLARK_HAVE_FFTWF || SKYLARK_HAVE_KISSFFT
    // Solve using Sylark
    timer.restart();
    skylark::nla::FasterLeastSquares(El::NORMAL, A, b, x, context);
    telp = timer.elapsed();
    check_solution(A, b, x, r, res, resAtr, resFac);
    if (rank == 0)
        std::cout << "Skylark:\t\t\t||r||_2 =  "
                  << boost::format("%.2f") % res
                  << " (x " << boost::format("%.5f") % (res / res_opt) << ")"
                  << "\t||r - r*||_2 / ||b - r*||_2 = " << boost::format("%.2e") % resFac
                  << "\t||A' * r||_2 = " << boost::format("%.2e") % resAtr
                  << "\t\tTime: " << boost::format("%.2e") % telp << " sec"
                  << std::endl;

    // Approximately solve using Sylark
    timer.restart();
    skylark::nla::ApproximateLeastSquares(El::NORMAL, A, b, x, context);
    telp = timer.elapsed();
    check_solution(A, b, x, r, res, resAtr, resFac);
    if (rank == 0)
        std::cout << "Skylark (approximate):\t\t||r||_2 =  "
                  << boost::format("%.2f") % res
                  << " (x " << boost::format("%.5f") % (res / res_opt) << ")"
                  << "\t||r - r*||_2 / ||b - r*||_2 = " << boost::format("%.2e") % resFac
                  << "\t||A' * r||_2 = " << boost::format("%.2e") % resAtr
                  << "\t\tTime: " << boost::format("%.2e") % telp << " sec"
                  << std::endl;
#else
    std::cout << "You need to have Skylark supporting FFTW or FFTWF " 
              << "to solve with skylark least_squares.cpp"
              << std::endl;
#endif
}



int main(int argc, char** argv) {

    El::Initialize(argc, argv);

    boost::mpi::communicator world;
    int rank = world.rank();

    if (rank == 0)
        std::cout << "Matrix: [VC,STAR], Rhs: [VC,STAR], Sol: [STAR,STAR]\n\n";
    experiment<El::DistMatrix<double, El::VC, El::STAR>,
               El::DistMatrix<double, El::VC, El::STAR>,
               El::DistMatrix<double, El::STAR, El::STAR> > ();

    if (rank == 0)
        std::cout << "\nMatrix: [MC,MR], Rhs: [MC,MR], Sol: [MC,MR]\n\n";
    experiment<El::DistMatrix<double>, El::DistMatrix<double>,
               El::DistMatrix<double> >();

    return 0;
}

For an example of sketch-and-solve regression in Python, run:

wget http://vikas.sindhwani.org/data.tar.gz
tar -xvzf data.tar.gz
mpiexec -np 4 python $SKYLARK_SRC_DIR/python-skylark/skylark/examples/sketch_solve.py data/usps.train data/usps.test

SVD

In $SKYLARK_SRC_DIR/nla/skylark_svd.cpp an example is provided that illustrates randomized singular value decompositions.

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
#include <El.hpp>

#include <boost/program_options.hpp>
#include <boost/mpi.hpp>
#include <boost/format.hpp>

#define SKYLARK_NO_ANY
#include <skylark.hpp>

#include <iostream>

namespace bpo = boost::program_options;

enum file_types {
    LIBSVM,
    ARC_LIST
};

template<typename InputType, typename FactorType, typename UType = FactorType,
         typename YType = FactorType>
void execute(bool directory, const std::string &fname,
    const std::string &hdfs, int port, const std::vector<int> &profile, int k,
    const skylark::nla::approximate_svd_params_t &params,
    const std::string &prefix,
    skylark::base::context_t &context) {

    boost::mpi::communicator world;
    int rank = world.rank();

    InputType A;
    FactorType S, V;
    UType U;
    YType Y;

    boost::mpi::timer timer;

    if (profile.empty()) {
        // Load A and Y (Y is thrown away)
        if (rank == 0) {
            std::cout << "Reading the matrix... ";
            std::cout.flush();
            timer.restart();
        }

        if (!hdfs.empty()) {
#       if SKYLARK_HAVE_LIBHDFS

            hdfsFS fs;
            if (rank == 0)
                fs = hdfsConnect(hdfs.c_str(), port);

            if (directory)
                SKYLARK_THROW_EXCEPTION(skylark::base::io_exception() <<
                    skylark::base::error_msg(
                                             "HDFS directory reading not yet supported."))
                else
                    skylark::utility::io::ReadLIBSVM(
                                                     fs, fname, A, Y, skylark::base::ROWS);

#       else

            SKYLARK_THROW_EXCEPTION(skylark::base::io_exception() <<
                skylark::base::error_msg("Install libhdfs for HDFS support!"));

#       endif
        } else {
            if (directory)
                skylark::utility::io::ReadDirLIBSVM(
                                                    fname, A, Y, skylark::base::ROWS);
            else
                skylark::utility::io::ReadLIBSVM(fname, A, Y, skylark::base::ROWS);
        }

        Y.Empty();
    } else {

        if (rank == 0) {
            std::cout << "Generating random matrix... ";
            std::cout.flush();
            timer.restart();
        }

        skylark::base::UniformMatrix(A, profile[0], profile[1], context);
    }

    if (rank == 0)
        std::cout <<"took " << boost::format("%.2e") % timer.elapsed()
                  << " sec\n";

    /* Compute approximate SVD */
    if (rank == 0) {
        std::cout << "Computing approximate SVD...";
        std::cout.flush();
        timer.restart();
    }

    skylark::nla::ApproximateSVD(A, U, S, V, k, context, params);

    if (rank == 0)
        std::cout <<"Took " << boost::format("%.2e") % timer.elapsed()
                  << " sec\n";

    /* Write results */
    if (rank == 0) {
        std::cout << "Writing results...";
        std::cout.flush();
        timer.restart();
    }

    El::Write(U, prefix + ".U", El::ASCII);
    El::Write(S, prefix + ".S", El::ASCII);
    El::Write(V, prefix + ".V", El::ASCII);

    if (rank == 0)
        std::cout <<"took " << boost::format("%.2e") % timer.elapsed()
                  << " sec\n";
}

template<typename InputType, typename FactorType, typename UType = FactorType,
         typename YType = FactorType>
void execute_sym(bool directory, const std::string &fname,
    const std::string& ftype,
    const std::string &hdfs, int port,
    const std::vector<int> &profile, bool lower, int k,
    const skylark::nla::approximate_svd_params_t &params,
    const std::string &prefix,
    skylark::base::context_t &context) {

    boost::mpi::communicator world;
    int rank = world.rank();

    InputType A;
    FactorType S, V;
    YType Y;

    boost::mpi::timer timer;

    if (profile.empty()) {
        // Load A and Y (Y is thrown away)
        if (rank == 0) {
            std::cout << "Reading the matrix... ";
            std::cout.flush();
            timer.restart();
        }

        if (!hdfs.empty()) {
#       if SKYLARK_HAVE_LIBHDFS

            hdfsFS fs;
            if (rank == 0)
                fs = hdfsConnect(hdfs.c_str(), port);

            if (directory)
                SKYLARK_THROW_EXCEPTION(skylark::base::io_exception() <<
                    skylark::base::error_msg(
                                             "HDFS directory reading not yet supported."))
                else
                    skylark::utility::io::ReadLIBSVM(
                                                     fs, fname, A, Y, skylark::base::ROWS);

#       else

            SKYLARK_THROW_EXCEPTION(skylark::base::io_exception() <<
                skylark::base::error_msg("Install libhdfs for HDFS support!"));

#       endif
        } else {
            // FIXME: ugly, fix options
            if (ftype.compare("ARC_LIST") == 0) {
                skylark::utility::io::ReadArcList(fname, A, world, true);
            } else {
                if (directory)
                    skylark::utility::io::ReadDirLIBSVM(
                                                        fname, A, Y, skylark::base::ROWS);
                else
                    skylark::utility::io::ReadLIBSVM(
                                                     fname, A, Y, skylark::base::ROWS);
            }
        }

        Y.Empty();
        
    } else {

        SKYLARK_THROW_EXCEPTION(skylark::base::unsupported_base_operation() <<
            skylark::base::error_msg("Uniform symmetric matrix generating not supported yet."));
    }

    if (rank == 0)
        std::cout << "took " << boost::format("%.2e") % timer.elapsed()
                  << " sec\n";

    /* Compute approximate SVD */
    if (rank == 0) {
        std::cout << "Computing approximate SVD...";
        std::cout.flush();
        timer.restart();
    }

    skylark::nla::ApproximateSymmetricSVD(lower ? El::LOWER : El::UPPER,
        A, V, S, k, context, params);

    if (rank == 0)
        std::cout <<"Took " << boost::format("%.2e") % timer.elapsed()
                  << " sec\n";

    /* Write results */
    if (rank == 0) {
        std::cout << "Writing results...";
        std::cout.flush();
        timer.restart();
    }

    El::Write(S, prefix + ".S", El::ASCII);
    El::Write(V, prefix + ".V", El::ASCII);

    if (rank == 0)
        std::cout <<"took " << boost::format("%.2e") % timer.elapsed()
                  << " sec\n";

    world.barrier();
}


int main(int argc, char* argv[]) {

    El::Initialize(argc, argv);

    boost::mpi::communicator world;
    int rank = world.rank();
    int size = world.size();

    int seed, k, powerits, port;
    std::string fname, ftype, prefix, hdfs;
    bool as_symmetric, as_sparse, skipqr, use_single, lower, directory;
    int oversampling_ratio, oversampling_additive;
    std::vector<int> profile;
    
    // Parse options
    bpo::options_description desc("Options");
    desc.add_options()
        ("help,h", "produce a help message")
        ("inputfile",
            bpo::value<std::string>(&fname),
            "Input file to run approximate SVD on (default libsvm format).")
        ("filetype",
            bpo::value<std::string>(&ftype),
            "Input file type (LIBSVM or ARC_LIST).")
        ("directory,d", "Whether inputfile is a directory of files whose"
            " concatination is the input.")
        ("seed,s",
            bpo::value<int>(&seed)->default_value(38734),
            "Seed for random number generation. OPTIONAL.")
        ("hdfs",
            bpo::value<std::string>(&hdfs)->default_value(""),
            "If not empty, will assume file is in an HDFS. "
            "Parameter is filesystem name.")
        ("port",
            bpo::value<int>(&port)->default_value(0),
            "For HDFS: port to use.")
        ("rank,k",
            bpo::value<int>(&k)->default_value(6),
            "Target rank. OPTIONAL.")
        ("powerits,i",
            bpo::value<int>(&powerits)->default_value(2),
            "Number of power iterations. OPTIONAL.")
        ("skipqr", "Whether to skip QR in each iteration. Higher than one power"
            " iterations is not recommended in this mode.")
        ("ratio,r",
            bpo::value<int>(&oversampling_ratio)->default_value(2),
            "Ratio of oversampling of rank. OPTIONAL.")
        ("additive,a",
            bpo::value<int>(&oversampling_additive)->default_value(0),
            "Additive factor for oversampling of rank. OPTIONAL.")
        ("symmetric", "Whether to treat the matrix as symmetric. "
         "Only upper part will be accessed, unless --lower is used.")
        ("lower", "For symmetric matrix, access only lower part (upper is default).")
        ("sparse", "Whether to load the matrix as a sparse one.")
        ("single", "Whether to use single precision instead of double.")
        ("profile",
            bpo::value<std::vector<int> >()->multitoken(),
            "Generate random matrix and run on it (for profiling)."
            "Requires specification of height and width. OPTIONAL.")
        ("prefix",
            bpo::value<std::string>(&prefix)->default_value("out"),
            "Prefix for output files (prefix.U.txt, prefix.S.txt"
            " and prefix.V.txt. OPTIONAL.");

    bpo::positional_options_description positional;
    positional.add("inputfile", 1);

    bpo::variables_map vm;
    try {
        bpo::store(bpo::command_line_parser(argc, argv)
            .options(desc).positional(positional).run(), vm);

        if (vm.count("help")) {
            if (rank == 0) {
                std::cout << "Usage: " << argv[0]
                          << " [options] input-file-name" << std::endl;
                std::cout << desc;
            }
            world.barrier();
            return 0;
        }

        if (vm["profile"].empty() && !vm.count("inputfile")) {
            if (rank == 0)
                std::cout << "Input file is required." << std::endl;
            world.barrier();
            return -1;
        }

        if (!vm["profile"].empty()) {
            if (vm["profile"].as<std::vector<int> >().size() < 2) {
                if (rank == 0)
                    std::cout << "Please specify height and width for --profile."
                              << std::endl;
                world.barrier();
                return -1;
            } else
                profile = vm["profile"].as<std::vector<int> >();
        }
        
        bpo::notify(vm);

        as_symmetric = vm.count("symmetric");
        as_sparse = vm.count("sparse");
        skipqr = vm.count("skipqr");
        use_single = vm.count("single");
        directory = vm.count("directory");
        lower = vm.count("lower");

    } catch(bpo::error& e) {
        if (rank == 0) {
            std::cerr << e.what() << std::endl;
            std::cerr << desc << std::endl;
        }
        world.barrier();
        return -1;
    }

    skylark::base::context_t context(seed);

    skylark::nla::approximate_svd_params_t params;
    params.skip_qr = skipqr;
    params.num_iterations = powerits;
    params.oversampling_ratio = oversampling_ratio;
    params.oversampling_additive = oversampling_additive;

    SKYLARK_BEGIN_TRY()

        if (size == 1) {
            if (!as_symmetric) {

                if (use_single) {
                    if (as_sparse)
                        execute<skylark::base::sparse_matrix_t<float>,
                                El::Matrix<float> >(
                                    directory, fname, hdfs,
                                    port, profile, k, params, prefix, context);
                    else
                        execute<El::Matrix<float>,
                                El::Matrix<float> >(directory, fname, hdfs,
                                    port, profile, k, params, prefix, context);

                } else {
                    if (as_sparse)
                        execute<skylark::base::sparse_matrix_t<double>,
                                El::Matrix<double> >(
                                     directory, fname, hdfs,
                                     port, profile, k, params, prefix, context);
                    else
                        execute<El::Matrix<double>,
                                El::Matrix<double> >(directory, fname, hdfs,
                                    port, profile, k, params, prefix, context);
                }

            } else {

                if (use_single) {
                    if (as_sparse)
                        execute_sym<skylark::base::sparse_matrix_t<float>,
                                    El::Matrix<float> >(
                                       directory, fname, ftype,
                                       hdfs, port, profile, lower, k, params, prefix,
                                        context);
                    else
                        execute_sym<El::Matrix<float>,
                                    El::Matrix<float> >(directory, fname, ftype,
                                        hdfs, port, profile, lower, k, params, prefix,
                                        context);

                } else {
                    if (as_sparse)
                        execute_sym<skylark::base::sparse_matrix_t<double>,
                                    El::Matrix<double>,
                                    El::Matrix<double> >(
                                        directory, fname, ftype, hdfs, port,
                                        profile, lower, k, params, prefix, context);
                    else
                        execute_sym<El::Matrix<double>,
                                    El::Matrix<double> >(directory, fname,
                                        ftype, hdfs, port, profile, lower, k, params,
                                        prefix, context);
                }
            }

        } else {
            if (!as_symmetric) {

                if (use_single) {
                    if (as_sparse)
                        execute<skylark::base::sparse_vc_star_matrix_t<float>,
                                El::DistMatrix<float, El::STAR, El::STAR>,
                                El::DistMatrix<float, El::VC, El::STAR>,
                                El::DistMatrix<float, El::VC, El::STAR> >(
                                    directory, fname, hdfs,
                                    port, profile, k, params, prefix, context);
                    else
                        execute<El::DistMatrix<float>,
                                El::DistMatrix<float> >(directory, fname, hdfs,
                                    port, profile, k, params, prefix, context);

                } else {
                    if (as_sparse)
                        execute<skylark::base::sparse_vc_star_matrix_t<double>,
                                El::DistMatrix<double, El::STAR, El::STAR>,
                                El::DistMatrix<double, El::VC, El::STAR>,
                                El::DistMatrix<double, El::VC, El::STAR> >(
                                    directory, fname, hdfs,
                                    port, profile, k, params, prefix, context);
                    else
                        execute<El::DistMatrix<double>,
                                El::DistMatrix<double> >(directory, fname, hdfs,
                                    port, profile, k, params, prefix, context);
                }

            } else {

                if (use_single) {
                    if (as_sparse)
                        execute_sym<skylark::base::sparse_vc_star_matrix_t<float>,
                                    El::DistMatrix<float, El::VC, El::STAR>,
                                    El::DistMatrix<float, El::VC, El::STAR> >(
                                        directory, fname, ftype,
                                        hdfs, port, profile, lower, k, params, prefix,
                                        context);
                    else
                        execute_sym<El::DistMatrix<float>,
                                    El::DistMatrix<float> >(directory, fname, ftype,
                                        hdfs, port, profile, lower, k, params, prefix,
                                        context);

                } else {
                    if (as_sparse)
                        execute_sym<skylark::base::sparse_vc_star_matrix_t<double>,
                                    El::DistMatrix<double, El::VC, El::STAR>,
                                    El::DistMatrix<double, El::VC, El::STAR> >(
                                        directory, fname, ftype, hdfs, port,
                                        profile, lower,
                                        k, params, prefix, context);
                    else
                        execute_sym<El::DistMatrix<double>,
                                    El::DistMatrix<double> >(directory, fname,
                                        ftype, hdfs, port, profile, lower, k, params,
                                        prefix, context);
                }
            }
        }


    SKYLARK_END_TRY() SKYLARK_CATCH_AND_PRINT((rank == 0))

    El::Finalize();
    return 0;
}

«  Introduction   ::   Contents   ::   Building from Source  »

© Copyright IBM Corporation 2012-2014. All Rights Reserved. Created using Sphinx 1.2.2.