diff --git a/skypy/galaxies/spectrum.py b/skypy/galaxies/spectrum.py
index 3d0722ed2..f54927610 100644
--- a/skypy/galaxies/spectrum.py
+++ b/skypy/galaxies/spectrum.py
@@ -143,6 +143,28 @@ def stellar_mass(self, coefficients, magnitudes, filter):
         Mt = self.absolute_magnitudes(coefficients, filter)
         return np.power(10, 0.4*(Mt-magnitudes))
 
+    def stellar_mass_remain(self, coefficients, magnitudes, filter):
+        r'''Compute total surviving stellar mass from absolute magnitudes.
+
+        Parameters
+        ----------
+        coefficients : (ng, nt) array_like
+            Array of template coefficients.
+        magnitudes : (ng,) array_like
+            The magnitudes to match in the reference bandpass.
+        filter : str
+            A single reference bandpass filter specification for
+            `~speclite.filters.load_filters`.
+
+        Returns
+        -------
+        stellar_mass : (ng,) array_like
+            Total surviving stellar mass of each galaxy in template units.
+        '''
+        m = self.stellar_mass(coefficients, magnitudes, filter)
+        m *= np.dot(coefficients, self.mremain)
+        return m
+
     def metallicity(self, coefficients):
         r'''Galaxy metallicities from kcorrect templates.
 
diff --git a/skypy/galaxies/tests/test_spectrum.py b/skypy/galaxies/tests/test_spectrum.py
index 5aea0f439..f52953060 100644
--- a/skypy/galaxies/tests/test_spectrum.py
+++ b/skypy/galaxies/tests/test_spectrum.py
@@ -158,6 +158,39 @@ def test_kcorrect_stellar_mass():
     np.testing.assert_allclose(stellar_mass, truth)
 
 
+@pytest.mark.skipif(not HAS_SPECLITE, reason='test requires speclite')
+def test_kcorrect_stellar_mass_remain():
+
+    from astropy import units
+    from skypy.galaxies.spectrum import kcorrect
+    from speclite.filters import FilterResponse
+
+    # Gaussian bandpass
+    filt_lam = np.logspace(3, 4, 1000) * units.AA
+    filt_mean = 5000 * units.AA
+    filt_width = 100 * units.AA
+    filt_tx = np.exp(-((filt_lam-filt_mean)/filt_width)**2)
+    filt_tx[[0, -1]] = 0
+    FilterResponse(wavelength=filt_lam, response=filt_tx,
+                   meta=dict(group_name='test', band_name='filt'))
+
+    # Using the identity matrix for the coefficients yields trivial test cases
+    coeff = np.eye(5)
+    Mt = kcorrect.absolute_magnitudes(coeff, 'test-filt')
+
+    # Using the absolute magnitudes of the templates as reference magnitudes
+    # should return one solar mass for each template.
+    stellar_mass = kcorrect.stellar_mass_remain(coeff, Mt, 'test-filt')
+    truth = kcorrect.mremain
+    np.testing.assert_allclose(stellar_mass, truth)
+
+    # Solution for given magnitudes without template mixing
+    Mb = np.array([10, 20, 30, 40, 50])
+    stellar_mass = kcorrect.stellar_mass_remain(coeff, Mb, 'test-filt')
+    truth = np.power(10, -0.4*(Mb-Mt)) * kcorrect.mremain
+    np.testing.assert_allclose(stellar_mass, truth)
+
+
 def test_kcorrect_metallicity():
 
     from skypy.galaxies.spectrum import kcorrect