diff --git a/brainiak/funcalign/rsrm.py b/brainiak/funcalign/rsrm.py
new file mode 100644
index 000000000..34efac464
--- /dev/null
+++ b/brainiak/funcalign/rsrm.py
@@ -0,0 +1,560 @@
+#  Copyright 2016 Intel Corporation
+#
+#  Licensed under the Apache License, Version 2.0 (the "License");
+#  you may not use this file except in compliance with the License.
+#  You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#  Unless required by applicable law or agreed to in writing, software
+#  distributed under the License is distributed on an "AS IS" BASIS,
+#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#  See the License for the specific language governing permissions and
+#  limitations under the License.
+"""Robust Shared Response Model (RSRM)
+
+The implementation is based on the following publications:
+
+.. [Turek2017] "Capturing Shared and Individual Information in fMRI Data",
+   J. Turek, C. Ellis, L. Skalaban, N. Turk-Browne, T. Willke
+   under review, 2017.
+"""
+
+# Authors: Javier Turek (Intel Labs), 2017
+
+import logging
+
+import numpy as np
+from sklearn.base import BaseEstimator, TransformerMixin
+from sklearn.utils import assert_all_finite
+from sklearn.utils.validation import NotFittedError
+
+__all__ = [
+    "RSRM"
+]
+
+logger = logging.getLogger(__name__)
+
+
+class RSRM(BaseEstimator, TransformerMixin):
+    """Robust Shared Response Model (RSRM)
+
+    Given multi-subject data, factorize it as a shared response R among all
+    subjects, an orthogonal transform W per subject, and an individual
+    (outlying) sparse component S per subject:
+
+    .. math:: X_i \\approx W_i R + S_i, \\forall i=1 \\dots N
+
+    This unsupervised model allows to learn idiosyncratic information for
+    subjects and simultaneously improve the shared response estimation.
+    The model has similar properties to the Shared Response Model (SRM) with
+    the addition of the individual components.
+
+    The model is estimated solving the following optimization problem:
+
+    .. math:: \\min_{W_i, S_i, R}\\sum_i \\frac{1}{2}\|X_i - W_i R - S_i\|_F^2
+    .. math:: + \\gamma\|S_i\|_1
+    .. math:: s.t. \\qquad W_i^TW_i = I \\quad \\forall i=1 \\dots N
+
+    The solution to this problem is obtained by applying a Block-Coordinate
+    Descent procedure. More details can be found in [Turek2017]_.
+
+    Parameters
+    ----------
+
+    n_iter : int, default: 10
+        Number of iterations to run the algorithm.
+
+    features : int, default: 50
+        Number of features to compute.
+
+    gamma : float, default: 1.0
+        Regularization parameter for the sparseness of the individual
+        components. Higher values yield sparser individual components.
+
+    rand_seed : int, default: 0
+        Seed for initializing the random number generator.
+
+
+    Attributes
+    ----------
+
+    w_ : list of array, element i has shape=[voxels_i, features]
+        The orthogonal transforms (mappings) for each subject.
+
+    r_ : array, shape=[features, timepoints]
+        The shared response.
+
+    s_ : list of array, element i has shape=[voxels_i, timepoints]
+        The individual components for each subject.
+
+    random_state_: `RandomState`
+        Random number generator initialized using rand_seed
+
+    Note
+    ----
+
+        The number of voxels may be different between subjects. However, the
+        number of timepoints for the alignment data must be the same across
+        subjects.
+
+        The Robust Shared Response Model is approximated using the
+        Block-Coordinate Descent (BCD) algorithm proposed in [Turek2017]_.
+
+        This is a single node version.
+    """
+
+    def __init__(self, n_iter=10, features=50, gamma=1.0, rand_seed=0):
+        self.n_iter = n_iter
+        self.features = features
+        self.lam = gamma
+        self.rand_seed = rand_seed
+
+    def fit(self, X):
+        """Compute the Robust Shared Response Model
+
+        Parameters
+        ----------
+
+        X : list of 2D arrays, element i has shape=[voxels_i, timepoints]
+            Each element in the list contains the fMRI data of one subject.
+        """
+        logger.info('Starting RSRM')
+
+        # Check that the regularizer value is positive
+        if 0.0 >= self.lam:
+            raise ValueError("Gamma parameter should be positive.")
+
+        # Check the number of subjects
+        if len(X) <= 1:
+            raise ValueError("There are not enough subjects in the input "
+                             "data to train the model.")
+
+        # Check for input data sizes
+        if X[0].shape[1] < self.features:
+            raise ValueError(
+                "There are not enough timepoints to train the model with "
+                "{0:d} features.".format(self.features))
+
+        # Check if all subjects have same number of TRs for alignment
+        number_trs = X[0].shape[1]
+        number_subjects = len(X)
+        for subject in range(number_subjects):
+            assert_all_finite(X[subject])
+            if X[subject].shape[1] != number_trs:
+                raise ValueError("Different number of alignment timepoints "
+                                 "between subjects.")
+
+        # Create a new random state
+        self.random_state_ = np.random.RandomState(self.rand_seed)
+
+        # Run RSRM
+        self.w_, self.r_, self.s_ = self._rsrm(X)
+
+        return self
+
+    def transform(self, X):
+        """Use the model to transform new data to Shared Response space
+
+        Parameters
+        ----------
+
+        X : list of 2D arrays, element i has shape=[voxels_i, timepoints_i]
+            Each element in the list contains the fMRI data of one subject.
+
+        Returns
+        -------
+
+        r : list of 2D arrays, element i has shape=[features_i, timepoints_i]
+            Shared responses from input data (X)
+
+        s : list of 2D arrays, element i has shape=[voxels_i, timepoints_i]
+            Individual data obtained from fitting model to input data (X)
+        """
+        # Check if the model exist
+        if hasattr(self, 'w_') is False:
+            raise NotFittedError("The model fit has not been run yet.")
+
+        # Check the number of subjects
+        if len(X) != len(self.w_):
+            raise ValueError("The number of subjects does not match the one"
+                             " in the model.")
+
+        r = [None] * len(X)
+        s = [None] * len(X)
+        for subject in range(len(X)):
+            if X[subject] is not None:
+                r[subject], s[subject] = self._transform_new_data(X[subject],
+                                                                  subject)
+
+        return r, s
+
+    def _transform_new_data(self, X, subject):
+        """Transform new data for a subjects by projecting to the shared subspace and
+        computing the individual information.
+
+        Parameters
+        ----------
+
+        X : array, shape=[voxels, timepoints]
+            The fMRI data of the subject.
+
+        subject : int
+            The subject id.
+
+        Returns
+        -------
+
+        R : array, shape=[features, timepoints]
+            Shared response from input data (X)
+
+        S : array, shape=[voxels, timepoints]
+            Individual data obtained from fitting model to input data (X)
+        """
+        S = np.zeros_like(X)
+        R = None
+        for i in range(self.n_iter):
+            R = self.w_[subject].T.dot(X - S)
+            S = self._shrink(X - self.w_[subject].dot(R), self.lam)
+        return R, S
+
+    def transform_subject(self, X):
+        """Transform a new subject using the existing model
+
+        Parameters
+        ----------
+
+        X : 2D array, shape=[voxels, timepoints]
+            The fMRI data of the new subject.
+
+        Returns
+        -------
+
+        w : 2D array, shape=[voxels, features]
+            Orthogonal mapping `W_{new}` for new subject
+
+        s : 2D array, shape=[voxels, timepoints]
+            Individual term `S_{new}` for new subject
+        """
+        # Check if the model exist
+        if hasattr(self, 'w_') is False:
+            raise NotFittedError("The model fit has not been run yet.")
+
+        # Check the number of TRs in the subject
+        if X.shape[1] != self.r_.shape[1]:
+            raise ValueError("The number of timepoints(TRs) does not match the"
+                             "one in the model.")
+
+        s = np.zeros_like(X)
+        for i in range(self.n_iter):
+            w = self._update_transform_subject(X, s, self.r_)
+            s = self._shrink(X - w.dot(self.r_), self.lam)
+
+        return w, s
+
+    def _rsrm(self, X):
+        """Block-Coordinate Descent algorithm for fitting RSRM.
+
+        Parameters
+        ----------
+
+        X : list of 2D arrays, element i has shape=[voxels_i, timepoints]
+            Each element in the list contains the fMRI data for alignment of
+            one subject.
+
+        Returns
+        -------
+
+        W : list of array, element i has shape=[voxels_i, features]
+            The orthogonal transforms (mappings) :math:`W_i` for each subject.
+
+        R : array, shape=[features, timepoints]
+            The shared response.
+
+        S : list of array, element i has shape=[voxels_i, timepoints]
+            The individual component :math:`S_i` for each subject.
+        """
+        subjs = len(X)
+        voxels = [X[i].shape[0] for i in range(subjs)]
+        TRs = X[0].shape[1]
+        features = self.features
+
+        # Initialization
+        W = self._init_transforms(subjs, voxels, features, self.random_state_)
+        S = self._init_individual(subjs, voxels, TRs)
+        R = self._update_shared_response(X, S, W, features)
+
+        if logger.isEnabledFor(logging.INFO):
+            objective = self._objective_function(X, W, R, S, self.lam)
+            logger.info('Objective function %f' % objective)
+
+        # Main loop
+        for i in range(self.n_iter):
+            W = self._update_transforms(X, S, R)
+            S = self._update_individual(X, W, R, self.lam)
+            R = self._update_shared_response(X, S, W, features)
+            # Print objective function every iteration
+            if logger.isEnabledFor(logging.INFO):
+                objective = self._objective_function(X, W, R, S, self.lam)
+                logger.info('Objective function %f' % objective)
+
+        return W, R, S
+
+    def _init_transforms(self, subjs, voxels, features, random_state):
+        """Initialize the mappings (Wi) with random orthogonal matrices.
+        Parameters
+        ----------
+
+        subjs : int
+            The number of subjects.
+
+        voxels : list of int
+            A list with the number of voxels per subject.
+
+        features : int
+            The number of features in the model.
+
+        random_state : `RandomState`
+            A random state to draw the mappings.
+
+        Returns
+        -------
+
+        W : list of array, element i has shape=[voxels_i, features]
+            The initialized orthogonal transforms (mappings) :math:`W_i` for
+            each subject.
+
+        Note
+        ----
+            Not thread safe.
+        """
+
+        # Init the Random seed generator
+        np.random.seed(self.rand_seed)
+
+        # Draw a random W for each subject
+        W = [random_state.random_sample((voxels[i], features))
+             for i in range(subjs)]
+        # Make it orthogonal it with QR decomposition
+        for i in range(subjs):
+            W[i], _ = np.linalg.qr(W[i])
+
+        return W
+
+    @staticmethod
+    def _objective_function(X, W, R, S, gamma):
+        """Evaluate the objective function.
+
+        .. math:: \sum_{i=1}^{N} 1/2 \| X_i - W_i R - S_i \|_F^2
+        .. math:: + /\gamma * \|S_i\|_1
+
+        Parameters
+        ----------
+
+        X : list of array, element i has shape=[voxels_i, timepoints]
+            Each element in the list contains the fMRI data for alignment of
+            one subject.
+
+        W : list of array, element i has shape=[voxels_i, features]
+            The orthogonal transforms (mappings) :math:`W_i` for each subject.
+
+        R : array, shape=[features, timepoints]
+            The shared response.
+
+        S : list of array, element i has shape=[voxels_i, timepoints]
+            The individual component :math:`S_i` for each subject.
+
+        gamma : float, default: 1.0
+            Regularization parameter for the sparseness of the individual
+            components.
+
+        Returns
+        -------
+
+        func : float
+            The RSRM objective function evaluated on the parameters to this
+            function.
+        """
+        subjs = len(X)
+        func = .0
+        for i in range(subjs):
+            func += 0.5 * np.sum((X[i] - W[i].dot(R) - S[i])**2) \
+                    + gamma * np.sum(np.abs(S[i]))
+        return func
+
+    @staticmethod
+    def _update_individual(X, W, R, gamma):
+        """Update the individual components `S_i`.
+
+        Parameters
+        ----------
+
+        X : list of 2D arrays, element i has shape=[voxels_i, timepoints]
+            Each element in the list contains the fMRI data for alignment of
+            one subject.
+
+        W : list of array, element i has shape=[voxels_i, features]
+            The orthogonal transforms (mappings) :math:`W_i` for each subject.
+
+        R : array, shape=[features, timepoints]
+            The shared response.
+
+        gamma : float, default: 1.0
+            Regularization parameter for the sparseness of the individual
+            components.
+
+        Returns
+        -------
+
+        S : list of array, element i has shape=[voxels_i, timepoints]
+            The individual component :math:`S_i` for each subject.
+        """
+        subjs = len(X)
+        S = []
+        for i in range(subjs):
+            S.append(RSRM._shrink(X[i] - W[i].dot(R), gamma))
+        return S
+
+    @staticmethod
+    def _init_individual(subjs, voxels, TRs):
+        """Initializes the individual components `S_i` to empty (all zeros).
+
+        Parameters
+        ----------
+
+        subjs : int
+            The number of subjects.
+
+        voxels : list of int
+            A list with the number of voxels per subject.
+
+        TRs : int
+            The number of timepoints in the data.
+
+        Returns
+        -------
+
+        S : list of 2D array, element i has shape=[voxels_i, timepoints]
+            The individual component :math:`S_i` for each subject initialized
+            to zero.
+        """
+        return [np.zeros((voxels[i], TRs)) for i in range(subjs)]
+
+    @staticmethod
+    def _update_shared_response(X, S, W, features):
+        """Update the shared response `R`.
+
+        Parameters
+        ----------
+
+        X : list of 2D arrays, element i has shape=[voxels_i, timepoints]
+            Each element in the list contains the fMRI data for alignment of
+            one subject.
+
+        S : list of array, element i has shape=[voxels_i, timepoints]
+            The individual component :math:`S_i` for each subject.
+
+        W : list of array, element i has shape=[voxels_i, features]
+            The orthogonal transforms (mappings) :math:`W_i` for each subject.
+
+        features : int
+            The number of features in the model.
+
+        Returns
+        -------
+
+        R : array, shape=[features, timepoints]
+            The updated shared response.
+        """
+        subjs = len(X)
+        TRs = X[0].shape[1]
+        R = np.zeros((features, TRs))
+        # Project the subject data with the individual component removed into
+        # the shared subspace and average over all subjects.
+        for i in range(subjs):
+            R += W[i].T.dot(X[i]-S[i])
+        R /= subjs
+        return R
+
+    @staticmethod
+    def _update_transform_subject(Xi, Si, R):
+        """Updates the mappings `W_i` for one subject.
+
+        Parameters
+        ----------
+
+        Xi : array, shape=[voxels, timepoints]
+            The fMRI data :math:`X_i` for aligning the subject.
+
+        Si : array, shape=[voxels, timepoints]
+            The individual component :math:`S_i` for the subject.
+
+        R : array, shape=[features, timepoints]
+            The shared response.
+
+        Returns
+        -------
+
+        Wi : array, shape=[voxels, features]
+            The orthogonal transform (mapping) :math:`W_i` for the subject.
+        """
+        A = Xi.dot(R.T)
+        A -= Si.dot(R.T)
+        # Solve the Procrustes problem
+        U, _, V = np.linalg.svd(A, full_matrices=False)
+        return U.dot(V)
+
+    @staticmethod
+    def _update_transforms(X, S, R):
+        """Updates the mappings `W_i` for each subject.
+
+        Parameters
+        ----------
+
+        X : list of 2D arrays, element i has shape=[voxels_i, timepoints]
+            Each element in the list contains the fMRI data for alignment of
+            one subject.ß
+
+        S : list of array, element i has shape=[voxels_i, timepoints]
+            The individual component :math:`S_i` for each subject.
+
+        R : array, shape=[features, timepoints]
+            The shared response.
+
+        Returns
+        -------
+
+        W : list of array, element i has shape=[voxels_i, features]
+            The orthogonal transforms (mappings) :math:`W_i` for each subject.
+        """
+        subjs = len(X)
+        W = []
+        for i in range(subjs):
+            W.append(RSRM._update_transform_subject(X[i], S[i], R))
+        return W
+
+    @staticmethod
+    def _shrink(v, gamma):
+        """Soft-shrinkage of an array with parameter gamma.
+
+        Parameters
+        ----------
+
+        v : array
+            Array containing the values to be applied to the shrinkage operator
+
+        gamma : float
+            Shrinkage parameter.
+
+        Returns
+        -------
+
+        v : array
+            The same input array after the shrinkage operator was applied.
+        """
+        pos = v > gamma
+        neg = v < -gamma
+        v[pos] -= gamma
+        v[neg] += gamma
+        v[np.logical_and(~pos, ~neg)] = .0
+        return v
diff --git a/docs/newsfragments/302.feature b/docs/newsfragments/302.feature
new file mode 100644
index 000000000..16c36c4f7
--- /dev/null
+++ b/docs/newsfragments/302.feature
@@ -0,0 +1 @@
+Added the Robust Shared Response Model method.
\ No newline at end of file
diff --git a/examples/funcalign/rsrm_synthetic_reconstruction.ipynb b/examples/funcalign/rsrm_synthetic_reconstruction.ipynb
new file mode 100644
index 000000000..b82ab4521
--- /dev/null
+++ b/examples/funcalign/rsrm_synthetic_reconstruction.ipynb
@@ -0,0 +1,191 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## This example replicates the synthetic experiment in the paper\n",
+    "## \"Capturing Shared and Individual Information in fMRI Data\", J. Turek, C. Ellis, L. Skalaban, N. Turk-Browne, T. Willke"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Import some libraries that we will need"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "import scipy.io\n",
+    "from scipy.stats import stats\n",
+    "import numpy as np\n",
+    "import brainiak.funcalign.srm as srm\n",
+    "import brainiak.funcalign.rsrm as rsrm\n",
+    "from mpl_toolkits.mplot3d import Axes3D\n",
+    "import matplotlib.pyplot as plt"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Experiment setup"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "voxels = 100\n",
+    "samples = 200\n",
+    "subjects = 10\n",
+    "features = 3\n",
+    "snr = 20 # in dB\n",
+    "amplitude = 8.0 # approximate amplitude used for the values in Si\n",
+    "p = 0.05 # probability of being non-zero in Si"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now we create some synthetic data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "# Create a Shared response R with K = 3\n",
+    "theta = np.linspace(-4 * np.pi, 4 * np.pi, samples)\n",
+    "z = np.linspace(-2, 2, samples)\n",
+    "r = z ** 2 + 1\n",
+    "x = r * np.sin(theta)\n",
+    "y = r * np.cos(theta)\n",
+    "curve = np.vstack((x, y, z))\n",
+    "print('Curve max, min values:', np.max(curve), np.min(curve))\n",
+    "\n",
+    "# Create the subjects' data\n",
+    "data = [None] * subjects\n",
+    "W = [None] * subjects\n",
+    "noise_level = 0.0\n",
+    "for s in range(subjects):\n",
+    "    R = curve\n",
+    "    W[s], _ = np.linalg.qr(np.random.randn(voxels, 3))\n",
+    "    data[s] = W[s].dot(R)\n",
+    "    noise_level += np.sum(np.abs(data[s])**2)\n",
+    "\n",
+    "# Compute noise_sigma from desired SNR\n",
+    "noise_level = noise_level / (10 ** (snr / 10))\n",
+    "noise_level = np.sqrt(noise_level  / subjects / voxels / samples)\n",
+    "\n",
+    "for s in range(subjects):\n",
+    "    n = noise_level * np.random.randn(voxels, samples)\n",
+    "    data[s] += n\n",
+    "    \n",
+    "S = [None] * subjects\n",
+    "for s in range(subjects):\n",
+    "    S[s] = (np.random.rand(data[s].shape[0], samples) < p) * ((np.random.rand(data[s].shape[0], samples) * amplitude - amplitude/2) )\n",
+    "    data[s] += S[s]"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Now we fit the algorithms, SRM and RSRM, to the synthetic data"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "algo_srm = srm.SRM(features=3, n_iter=20)\n",
+    "algo_srm.fit(data)\n",
+    "\n",
+    "\n",
+    "algo_rsrm = rsrm.RSRM(features=3, gamma=0.35, n_iter=20)\n",
+    "algo_rsrm.fit(data)\n",
+    "\n",
+    "print('Done')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The following function finds the orthogonal transform to align the shared response to the original curve."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "def find_orthogonal_transform(shared_response, curve):\n",
+    "    u,_,vt = np.linalg.svd(shared_response.dot(curve.T))\n",
+    "    q = u.dot(vt)\n",
+    "    aligned_curve = q.T.dot(shared_response)\n",
+    "    return aligned_curve"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Plot the results"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "fig = plt.figure()\n",
+    "ax = fig.gca(projection='3d')\n",
+    "ax.plot(curve[0, :], curve[1, :], curve[2, :], '-g', label='original', lineWidth=5)\n",
+    "\n",
+    "proj = find_orthogonal_transform(algo_srm.s_, curve)\n",
+    "ax.plot(proj[0, :], proj[1, :], proj[2, :], '-b', label='SRM', lineWidth=3)\n",
+    "proj = find_orthogonal_transform(algo_rsrm.r_, curve)\n",
+    "ax.plot(proj[0, :], proj[1, :], proj[2, :], '-r', label='RSRM', lineWidth=3)\n",
+    "plt.legend()"
+   ]
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "codemirror_mode": {
+    "name": "ipython",
+    "version": 3
+   },
+   "file_extension": ".py",
+   "mimetype": "text/x-python",
+   "name": "python",
+   "nbconvert_exporter": "python",
+   "pygments_lexer": "ipython3",
+   "version": "3.6.3"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 1
+}
diff --git a/examples/funcalign/srm_image_prediction_example.ipynb b/examples/funcalign/srm_image_prediction_example.ipynb
index 52610f3d9..378a50173 100644
--- a/examples/funcalign/srm_image_prediction_example.ipynb
+++ b/examples/funcalign/srm_image_prediction_example.ipynb
@@ -404,7 +404,7 @@
     "    cm[subject] = confusion_matrix(test_labels, predicted_labels)\n",
     "\n",
     "    # Normalize the confusion matrix\n",
-    "    cm[subject] = cm[subject].astype('float') / cm[subject].sum(axis=1)[:, np.newaxis]\n"
+    "    cm[subject] = cm[subject].astype('float') / cm[subject].sum(axis=1)[:, np.newaxis]"
    ]
   },
   {
@@ -436,7 +436,7 @@
   "language_info": {
    "codemirror_mode": {
     "name": "ipython",
-    "version": 3
+    "version": 3.0
    },
    "file_extension": ".py",
    "mimetype": "text/x-python",
@@ -448,4 +448,4 @@
  },
  "nbformat": 4,
  "nbformat_minor": 0
-}
+}
\ No newline at end of file
diff --git a/tests/funcalign/test_rsrm.py b/tests/funcalign/test_rsrm.py
new file mode 100644
index 000000000..b723a58c7
--- /dev/null
+++ b/tests/funcalign/test_rsrm.py
@@ -0,0 +1,130 @@
+#  Copyright 2016 Intel Corporation
+#
+#  Licensed under the Apache License, Version 2.0 (the "License");
+#  you may not use this file except in compliance with the License.
+#  You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#  Unless required by applicable law or agreed to in writing, software
+#  distributed under the License is distributed on an "AS IS" BASIS,
+#  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#  See the License for the specific language governing permissions and
+#  limitations under the License.
+from sklearn.exceptions import NotFittedError
+import pytest
+
+
+def test_can_instantiate():
+    import brainiak.funcalign.rsrm
+    s = brainiak.funcalign.rsrm.RSRM()
+    assert s, "Invalid RSRM instance!"
+
+    import numpy as np
+    np.random.seed(0)
+
+    voxels = 100
+    samples = 500
+    subjects = 2
+    features = 3
+
+    s = brainiak.funcalign.rsrm.RSRM(n_iter=5, features=features)
+    assert s, "Invalid RSRM instance!"
+
+    # Create a Shared response S with K = 3
+    theta = np.linspace(-4 * np.pi, 4 * np.pi, samples)
+    z = np.linspace(-2, 2, samples)
+    r = z**2 + 1
+    x = r * np.sin(theta)
+    y = r * np.cos(theta)
+
+    S = np.vstack((x, y, z))
+
+    X = []
+    W = []
+    Q, R = np.linalg.qr(np.random.random((voxels, features)))
+    W.append(Q)
+    X.append(Q.dot(S) + 0.1*np.random.random((voxels, samples)))
+
+    # Check that transform does NOT run before fitting the model
+    with pytest.raises(NotFittedError):
+        s.transform(X)
+    print("Test: transforming before fitting the model")
+
+    # Check that it does NOT run with 1 subject
+    with pytest.raises(ValueError):
+        s.fit(X)
+    print("Test: running RSRM with 1 subject")
+
+    for subject in range(1, subjects):
+        Q, R = np.linalg.qr(np.random.random((voxels, features)))
+        W.append(Q)
+        X.append(Q.dot(S) + 0.1*np.random.random((voxels, samples)))
+
+    # Check that runs with 2 subject
+    s.fit(X)
+
+    assert len(s.w_) == subjects, (
+        "Invalid computation of RSRM! (wrong # subjects in W)")
+    for subject in range(subjects):
+        assert s.s_[subject].shape[0] == voxels, (
+            "Invalid computation of RSRM! (wrong # voxels in S)")
+        assert s.s_[subject].shape[1] == samples, (
+            "Invalid computation of RSRM! (wrong # samples in S)")
+        assert s.w_[subject].shape[0] == voxels, (
+            "Invalid computation of RSRM! (wrong # voxels in W)")
+        assert s.w_[subject].shape[1] == features, (
+            "Invalid computation of RSRM! (wrong # features in W)")
+        ortho = np.linalg.norm(s.w_[subject].T.dot(s.w_[subject])
+                               - np.eye(s.w_[subject].shape[1]),
+                               'fro')
+        assert ortho < 1e-7, "A Wi mapping is not orthonormal in RSRM."
+        difference = np.linalg.norm(X[subject] - s.w_[subject].dot(s.r_),
+                                    'fro')
+        datanorm = np.linalg.norm(X[subject], 'fro')
+        assert difference/datanorm < 1.0, "Model seems incorrectly computed."
+    assert s.r_.shape[0] == features, (
+        "Invalid computation of RSRM! (wrong # features in R)")
+    assert s.r_.shape[1] == samples, (
+        "Invalid computation of RSRM! (wrong # samples in R)")
+
+    # Check that it does run to compute the shared response after the model
+    # computation
+    new_r, _ = s.transform(X)
+
+    assert len(new_r) == subjects, (
+        "Invalid computation of RSRM! (wrong # subjects after transform)")
+    for subject in range(subjects):
+        assert new_r[subject].shape[0] == features, (
+            "Invalid computation of RSRM! (wrong # features after transform)")
+        assert new_r[subject].shape[1] == samples, (
+            "Invalid computation of RSRM! (wrong # samples after transform)")
+
+    # Check that it does run to compute a new subject
+    new_w, new_s = s.transform_subject(X[0])
+    assert new_w.shape[1] == features, (
+            "Invalid computation of RSRM! (wrong # features for new subject)")
+    assert new_s.shape[1] == samples, (
+            "Invalid computation of RSRM! (wrong # samples for new subject)")
+    assert new_s.shape[0] == voxels, (
+            "Invalid computation of RSRM! (wrong # voxels for new subject)")
+    assert new_w.shape[0] == voxels, (
+            "Invalid computation of RSRM! (wrong # voxels for new subject)")
+
+    # Check that it does NOT run with non-matching number of subjects
+    with pytest.raises(ValueError):
+        s.transform(X[1])
+    print("Test: transforming with non-matching number of subjects")
+
+    # Check that it does not run without enough samples (TRs).
+    with pytest.raises(ValueError):
+        s.set_params(features=(samples+1))
+        s.fit(X)
+    print("Test: not enough samples")
+
+    # Check that it does not run with different number of samples (TRs)
+    S2 = S[:, :-2]
+    X.append(Q.dot(S2))
+    with pytest.raises(ValueError):
+        s.fit(X)
+    print("Test: different number of samples per subject")