diff --git a/nimfa/methods/factorization/bd.py b/nimfa/methods/factorization/bd.py
index b8a1e61..f3caea1 100644
--- a/nimfa/methods/factorization/bd.py
+++ b/nimfa/methods/factorization/bd.py
@@ -187,10 +187,10 @@ def __init__(self, V, seed=None, W=None, H=None, rank=30, max_iter=30,
         nmf_std.Nmf_std.__init__(self, vars())
         if self.alpha is None:
             self.alpha = sp.csr_matrix((self.V.shape[0], self.rank))
-        self.alpha = self.alpha.tocsr() if sp.isspmatrix(self.alpha) else np.mat(self.alpha)
+        self.alpha = self.alpha.tocsr() if sp.isspmatrix(self.alpha) else np.asmatrix(self.alpha)
         if self.beta is None:
             self.beta = sp.csr_matrix((self.rank, self.V.shape[1]))
-        self.beta = self.beta.tocsr() if sp.isspmatrix(self.beta) else np.mat(self.beta)
+        self.beta = self.beta.tocsr() if sp.isspmatrix(self.beta) else np.asmatrix(self.beta)
         if self.n_w is None:
             self.n_w = np.zeros((self.rank, 1))
         if self.n_h is None:
diff --git a/nimfa/methods/factorization/icm.py b/nimfa/methods/factorization/icm.py
index 89e6ab7..8abf6e4 100644
--- a/nimfa/methods/factorization/icm.py
+++ b/nimfa/methods/factorization/icm.py
@@ -162,10 +162,10 @@ def __init__(self, V, seed=None, W=None, H=None, H1=None,
         nmf_std.Nmf_std.__init__(self, vars())
         if self.alpha is None:
             self.alpha = sp.rand(self.V.shape[0], self.rank, density=0.8, format='csr')
-        self.alpha= self.alpha.tocsr() if sp.isspmatrix(self.alpha) else np.mat(self.alpha)
+        self.alpha= self.alpha.tocsr() if sp.isspmatrix(self.alpha) else np.asmatrix(self.alpha)
         if self.beta is None:
             self.beta = sp.rand(self.rank, self.V.shape[1], density=0.8, format='csr')
-        self.beta = self.beta.tocsr() if sp.isspmatrix(self.beta) else np.mat(self.beta)
+        self.beta = self.beta.tocsr() if sp.isspmatrix(self.beta) else np.asmatrix(self.beta)
         self.tracker = mf_track.Mf_track() if self.track_factor and self.n_run > 1 \
                                               or self.track_error else None
 
diff --git a/nimfa/methods/factorization/lfnmf.py b/nimfa/methods/factorization/lfnmf.py
index fb4a2a7..4530118 100644
--- a/nimfa/methods/factorization/lfnmf.py
+++ b/nimfa/methods/factorization/lfnmf.py
@@ -157,8 +157,8 @@ def factorize(self):
         for run in range(self.n_run):
             self.W, self.H = self.seed.initialize(
                 self.V, self.rank, self.options)
-            self.Sw, self.Sb = np.mat(
-                np.zeros((1, 1))), np.mat(np.zeros((1, 1)))
+            self.Sw, self.Sb = np.asmatrix(
+                np.zeros((1, 1))), np.asmatrix(np.zeros((1, 1)))
             p_obj = c_obj = sys.float_info.max
             best_obj = c_obj if run == 0 else best_obj
             iter = 0
@@ -244,7 +244,7 @@ def update(self):
         # update within class scatter and between class
         self.Sw = sum(sum(dot(self.H[:, c2m[i][j]] - avgs[i], (self.H[:, c2m[i][j]] - avgs[i]).T)
                           for j in range(len(c2m[i]))) for i in c2m)
-        avgs_t = np.mat(np.zeros((self.rank, 1)))
+        avgs_t = np.asmatrix(np.zeros((self.rank, 1)))
         for k in avgs:
             avgs_t += avgs[k]
         avgs_t /= len(avgs)
@@ -259,7 +259,7 @@ def _encoding(self, idxH):
             c2m.setdefault(idxH[0, i], [])
             c2m[idxH[0, i]].append(i)
             # compute mean value of class idx in encoding matrix H
-            avgs.setdefault(idxH[0, i], np.mat(np.zeros((self.rank, 1))))
+            avgs.setdefault(idxH[0, i], np.asmatrix(np.zeros((self.rank, 1))))
             avgs[idxH[0, i]] += self.H[:, i]
         for k in avgs:
             avgs[k] /= len(c2m[k])
diff --git a/nimfa/methods/factorization/nmf.py b/nimfa/methods/factorization/nmf.py
index 8c7d6cb..ca1ec98 100644
--- a/nimfa/methods/factorization/nmf.py
+++ b/nimfa/methods/factorization/nmf.py
@@ -303,7 +303,7 @@ def conn(self):
         cons = elop(mat1, mat2, eq)
         if not hasattr(self, 'consold'):
             self.cons = cons
-            self.consold = np.mat(np.logical_not(cons))
+            self.consold = np.asmatrix(np.logical_not(cons))
         else:
             self.consold = self.cons
             self.cons = cons
diff --git a/nimfa/methods/factorization/pmfcc.py b/nimfa/methods/factorization/pmfcc.py
index c5a7a07..e8aae10 100644
--- a/nimfa/methods/factorization/pmfcc.py
+++ b/nimfa/methods/factorization/pmfcc.py
@@ -133,7 +133,7 @@ def __init__(self, V, seed=None, W=None, H=None, rank=30, max_iter=30,
         self.aseeds = ["random", "fixed", "nndsvd", "random_c", "random_vcol"]
         smf.Smf.__init__(self, vars())
         if not self.Theta:
-            self.Theta = np.mat(np.zeros((self.V.shape[1], self.V.shape[1])))
+            self.Theta = np.asmatrix(np.zeros((self.V.shape[1], self.V.shape[1])))
         self.tracker = mf_track.Mf_track() if self.track_factor and self.n_run > 1 \
                                               or self.track_error else None
 
diff --git a/nimfa/methods/factorization/psmf.py b/nimfa/methods/factorization/psmf.py
index ef872eb..1fe7129 100644
--- a/nimfa/methods/factorization/psmf.py
+++ b/nimfa/methods/factorization/psmf.py
@@ -188,8 +188,8 @@ def factorize(self):
                 self.H = self.V.__class__(
                     (self.rank, self.V.shape[1]), dtype='d')
             else:
-                self.W = np.mat(np.zeros((self.V.shape[0], self.rank)))
-                self.H = np.mat(np.zeros((self.rank, self.V.shape[1])))
+                self.W = np.asmatrix(np.zeros((self.V.shape[0], self.rank)))
+                self.H = np.asmatrix(np.zeros((self.rank, self.V.shape[1])))
             self.s = np.zeros((self.V.shape[0], self.N), int)
             self.r = np.zeros((self.V.shape[0], 1), int)
             self.psi = np.array(std(self.V, axis=1, ddof=0))
diff --git a/nimfa/methods/factorization/sepnmf.py b/nimfa/methods/factorization/sepnmf.py
index 4b8db0e..7274e91 100644
--- a/nimfa/methods/factorization/sepnmf.py
+++ b/nimfa/methods/factorization/sepnmf.py
@@ -318,7 +318,7 @@ def norm_axis(X, p=None, axis=None):
     if sp.isspmatrix(X):
         return sla.norm(X, ord=p, axis=axis)
     else:
-        n = nla.norm(np.mat(X), ord=p, axis=axis)
+        n = nla.norm(np.asmatrix(X), ord=p, axis=axis)
         if axis is None:
             return n
         else:
diff --git a/nimfa/methods/factorization/snmf.py b/nimfa/methods/factorization/snmf.py
index c30d29b..215a8f3 100644
--- a/nimfa/methods/factorization/snmf.py
+++ b/nimfa/methods/factorization/snmf.py
@@ -183,8 +183,8 @@ def factorize(self):
             if sp.isspmatrix(self.H):
                 self.H = self.H.tolil()
             iter = 0
-            self.idx_w_old = np.mat(np.zeros((self.V.shape[0], 1)))
-            self.idx_h_old = np.mat(np.zeros((1, self.V.shape[1])))
+            self.idx_w_old = np.asmatrix(np.zeros((self.V.shape[0], 1)))
+            self.idx_h_old = np.asmatrix(np.zeros((1, self.V.shape[1])))
             c_obj = sys.float_info.max
             best_obj = c_obj if run == 0 else best_obj
             # count the number of convergence checks that column clusters and
@@ -200,7 +200,7 @@ def factorize(self):
                     sp.eye(self.rank, self.rank, format='lil')
             else:
                 self.beta_vec = sqrt(self.beta) * np.ones((1, self.rank))
-                self.I_k = self.eta * np.mat(np.eye(self.rank))
+                self.I_k = self.eta * np.asmatrix(np.eye(self.rank))
             self.n_restart = 0
             if self.callback_init:
                 self.final_obj = c_obj
@@ -287,8 +287,8 @@ def update(self):
             if self.n_restart >= 100:
                 raise utils.MFError(
                     "Too many restarts due to too large beta parameter.")
-            self.idx_w_old = np.mat(np.zeros((self.V.shape[0], 1)))
-            self.idx_h_old = np.mat(np.zeros((1, self.V.shape[1])))
+            self.idx_w_old = np.asmatrix(np.zeros((self.V.shape[0], 1)))
+            self.idx_h_old = np.asmatrix(np.zeros((1, self.V.shape[1])))
             self.inc = 0
             self.W, _ = self.seed.initialize(self.V, self.rank, self.options)
             # normalize W and convert to lil
@@ -333,7 +333,7 @@ def objective(self):
         resmat1 = elop(self.W, WHHt - VHt + self.eta ** 2 * self.W, min)
         res_vec = nz_data(resmat) + nz_data(resmat1)
         # L1 norm
-        self.conv = norm(np.mat(res_vec), 1)
+        self.conv = norm(np.asmatrix(res_vec), 1)
         err_avg = self.conv / len(res_vec)
         self.idx_w_old = idx_w
         self.idx_h_old = idx_h
@@ -388,7 +388,7 @@ def _spfcnnls(self, C, A):
             # make infeasible solutions feasible (standard NNLS inner loop)
             if len(h_set) > 0:
                 n_h_set = len(h_set)
-                alpha = np.mat(np.zeros((l_var, n_h_set)))
+                alpha = np.asmatrix(np.zeros((l_var, n_h_set)))
                 while len(h_set) > 0 and iter < max_iter:
                     iter += 1
                     alpha[:, :n_h_set] = np.Inf
@@ -465,15 +465,15 @@ def __spcssls(self, CtC, CtA, p_set=None):
             # equivalent if CtC is square matrix
             for k in range(CtA.shape[1]):
                 ls = sp.linalg.gmres(CtC, CtA[:, k].toarray())[0]
-                K[:, k] = sp.lil_matrix(np.mat(ls).T)
+                K[:, k] = sp.lil_matrix(np.asmatrix(ls).T)
             # K = dot(np.linalg.pinv(CtC), CtA)
         else:
             l_var, p_rhs = p_set.shape
             coded_p_set = dot(
-                sp.lil_matrix(np.mat(2 ** np.array(list(range(l_var - 1, -1, -1))))), p_set)
+                sp.lil_matrix(np.asmatrix(2 ** np.array(list(range(l_var - 1, -1, -1))))), p_set)
             sorted_p_set, sorted_idx_set = sort(coded_p_set.todense())
-            breaks = diff(np.mat(sorted_p_set))
-            break_idx = [-1] + find(np.mat(breaks)) + [p_rhs]
+            breaks = diff(np.asmatrix(sorted_p_set))
+            break_idx = [-1] + find(np.asmatrix(breaks)) + [p_rhs]
             for k in range(len(break_idx) - 1):
                 cols2solve = sorted_idx_set[
                     break_idx[k] + 1: break_idx[k + 1] + 1]
@@ -483,7 +483,7 @@ def __spcssls(self, CtC, CtA, p_set=None):
                 sol = sp.lil_matrix(K.shape)
                 for k in range(tmp_ls.shape[1]):
                     ls = sp.linalg.gmres(CtC[:, vars][vars, :], tmp_ls[:, k].toarray())[0]
-                    sol[:, k] = sp.lil_matrix(np.mat(ls).T)
+                    sol[:, k] = sp.lil_matrix(np.asmatrix(ls).T)
                 i = 0
                 for c in cols2solve:
                     j = 0
@@ -519,7 +519,7 @@ def _fcnnls(self, C, A):
         A = A.todense() if sp.isspmatrix(A) else A
         _, l_var = C.shape
         p_rhs = A.shape[1]
-        W = np.mat(np.zeros((l_var, p_rhs)))
+        W = np.asmatrix(np.zeros((l_var, p_rhs)))
         iter = 0
         max_iter = 3 * l_var
         # precompute parts of pseudoinverse
@@ -542,7 +542,7 @@ def _fcnnls(self, C, A):
             # make infeasible solutions feasible (standard NNLS inner loop)
             if len(h_set) > 0:
                 n_h_set = len(h_set)
-                alpha = np.mat(np.zeros((l_var, n_h_set)))
+                alpha = np.asmatrix(np.zeros((l_var, n_h_set)))
                 while len(h_set) > 0 and iter < max_iter:
                     iter += 1
                     alpha[:, :n_h_set] = np.Inf
@@ -597,7 +597,7 @@ def __cssls(self, CtC, CtA, p_set=None):
         Solve the set of equations CtA = CtC * K for variables defined in set p_set
         using the fast combinatorial approach (van Benthem and Keenan, 2004).
         """
-        K = np.mat(np.zeros(CtA.shape))
+        K = np.asmatrix(np.zeros(CtA.shape))
         if p_set is None or p_set.size == 0 or all(p_set):
             # equivalent if CtC is square matrix
             K = np.linalg.lstsq(CtC, CtA, rcond=-1)[0]
@@ -605,10 +605,10 @@ def __cssls(self, CtC, CtA, p_set=None):
         else:
             l_var, p_rhs = p_set.shape
             coded_p_set = dot(
-                np.mat(2 ** np.array(list(range(l_var - 1, -1, -1)))), p_set)
+                np.asmatrix(2 ** np.array(list(range(l_var - 1, -1, -1)))), p_set)
             sorted_p_set, sorted_idx_set = sort(coded_p_set)
-            breaks = diff(np.mat(sorted_p_set))
-            break_idx = [-1] + find(np.mat(breaks)) + [p_rhs]
+            breaks = diff(np.asmatrix(sorted_p_set))
+            break_idx = [-1] + find(np.asmatrix(breaks)) + [p_rhs]
             for k in range(len(break_idx) - 1):
                 cols2solve = sorted_idx_set[
                     break_idx[k] + 1: break_idx[k + 1] + 1]
diff --git a/nimfa/methods/factorization/snmnmf.py b/nimfa/methods/factorization/snmnmf.py
index c528511..e314c6b 100644
--- a/nimfa/methods/factorization/snmnmf.py
+++ b/nimfa/methods/factorization/snmnmf.py
@@ -192,10 +192,10 @@ def __init__(self, V, V1, seed=None, W=None, H=None, H1=None,
         nmf_mm.Nmf_mm.__init__(self, vars())
         if self.A is None:
             self.A = sp.csr_matrix((self.V1.shape[1], self.V1.shape[1]))
-        self.A = self.A.tocsr() if sp.isspmatrix(self.A) else np.mat(self.A)
+        self.A = self.A.tocsr() if sp.isspmatrix(self.A) else np.asmatrix(self.A)
         if self.B is None:
             self.B = sp.csr_matrix((self.V.shape[1], self.V1.shape[1]))
-        self.B = self.B.tocsr() if sp.isspmatrix(self.B) else np.mat(self.B)
+        self.B = self.B.tocsr() if sp.isspmatrix(self.B) else np.asmatrix(self.B)
         self.tracker = mf_track.Mf_track() if self.track_factor and self.n_run > 1 \
                                               or self.track_error else None