diff --git a/.gitignore b/.gitignore
index a30b5eb4..620482ae 100644
--- a/.gitignore
+++ b/.gitignore
@@ -31,6 +31,14 @@ coverage.xml
 ehthumbs.db
 Thumbs.db
 
+# virtual environments #
+########################
+.env
+.venv/
+.conda/
+.mamba/
+.python-version
+
 # Directories
 #############
 dist
diff --git a/pygam/pygam.py b/pygam/pygam.py
index bebcef13..a5c92293 100644
--- a/pygam/pygam.py
+++ b/pygam/pygam.py
@@ -703,7 +703,7 @@ def _initial_estimate(self, y, modelmat):
 
         # solve the linear problem
         return np.linalg.solve(
-            load_diagonal(modelmat.T.dot(modelmat).A), modelmat.T.dot(y_)
+            load_diagonal(modelmat.T.dot(modelmat).toarray()), modelmat.T.dot(y_)
         )
 
         # not sure if this is faster...
@@ -780,7 +780,7 @@ def _pirls(self, X, Y, weights):
             self._on_loop_start(vars())
 
             WB = W.dot(modelmat[mask, :])  # common matrix product
-            Q, R = np.linalg.qr(WB.A)
+            Q, R = np.linalg.qr(WB.toarray())
 
             if not np.isfinite(Q).all() or not np.isfinite(R).all():
                 raise ValueError(
@@ -1401,7 +1401,7 @@ def _get_quantiles(
         idxs = self.terms.get_coef_indices(term)
         cov = self.statistics_['cov'][idxs][:, idxs]
 
-        var = (modelmat.dot(cov) * modelmat.A).sum(axis=1)
+        var = (modelmat.dot(cov) * modelmat.toarray()).sum(axis=1)
         if prediction:
             var += self.distribution.scale
 
diff --git a/pygam/terms.py b/pygam/terms.py
index c033e726..83fb22af 100644
--- a/pygam/terms.py
+++ b/pygam/terms.py
@@ -1500,7 +1500,7 @@ def _build_marginal_constraints(self, i, coef, constraint_lam, constraint_l2):
             )
 
             # now enter it into the composite
-            composite_C[tuple(np.meshgrid(slice_, slice_))] = slice_C.A
+            composite_C[tuple(np.meshgrid(slice_, slice_))] = slice_C.toarray()
 
         return sp.sparse.csc_matrix(composite_C)
 
diff --git a/pygam/tests/test_penalties.py b/pygam/tests/test_penalties.py
index bf6dd68a..cb0a4193 100644
--- a/pygam/tests/test_penalties.py
+++ b/pygam/tests/test_penalties.py
@@ -23,13 +23,13 @@ def test_single_spline_penalty():
     monotonic_ and convexity_ should be 0.
     """
     coef = np.array(1.0)
-    assert np.alltrue(derivative(1, coef).A == 0.0)
-    assert np.alltrue(l2(1, coef).A == 1.0)
-    assert np.alltrue(monotonic_inc(1, coef).A == 0.0)
-    assert np.alltrue(monotonic_dec(1, coef).A == 0.0)
-    assert np.alltrue(convex(1, coef).A == 0.0)
-    assert np.alltrue(concave(1, coef).A == 0.0)
-    assert np.alltrue(none(1, coef).A == 0.0)
+    assert np.all(derivative(1, coef).toarray() == 0.0)
+    assert np.all(l2(1, coef).toarray() == 1.0)
+    assert np.all(monotonic_inc(1, coef).toarray() == 0.0)
+    assert np.all(monotonic_dec(1, coef).toarray() == 0.0)
+    assert np.all(convex(1, coef).toarray() == 0.0)
+    assert np.all(concave(1, coef).toarray() == 0.0)
+    assert np.all(none(1, coef).toarray() == 0.0)
 
 
 def test_wrap_penalty():
@@ -43,12 +43,12 @@ def test_wrap_penalty():
 
     fit_linear = True
     p = wrap_penalty(none, fit_linear, linear_penalty=linear_penalty)
-    P = p(n, coef).A
+    P = p(n, coef).toarray()
     assert P.sum() == linear_penalty
 
     fit_linear = False
     p = wrap_penalty(none, fit_linear, linear_penalty=linear_penalty)
-    P = p(n, coef).A
+    P = p(n, coef).toarray()
     assert P.sum() == 0.0
 
 
diff --git a/pygam/tests/test_terms.py b/pygam/tests/test_terms.py
index 72c9ce0b..a773ac99 100644
--- a/pygam/tests/test_terms.py
+++ b/pygam/tests/test_terms.py
@@ -315,10 +315,10 @@ def der1(*args, **kwargs):
 
     # check all the dimensions
     for i in range(3):
-        P = term._build_marginal_penalties(i).A
+        P = term._build_marginal_penalties(i).toarray()
         C = term._build_marginal_constraints(
             i, -np.arange(term.n_coefs), constraint_lam=1, constraint_l2=0
-        ).A
+        ).toarray()
 
         assert (P == C).all()
 
@@ -362,11 +362,11 @@ def test_compose_penalties(self):
         term = SplineTerm(feature=0, penalties=['auto', 'none'])
 
         # penalties should be equivalent
-        assert (term.build_penalties() == base_term.build_penalties()).A.all()
+        assert (term.build_penalties() == base_term.build_penalties()).toarray().all()
 
         # multitple penalties should be additive, not multiplicative,
         # so 'none' penalty should have no effect
-        assert np.abs(term.build_penalties().A).sum() > 0
+        assert np.abs(term.build_penalties().toarray()).sum() > 0
 
     def test_compose_constraints(self, hepatitis_X_y):
         """we should be able to compose penalties
diff --git a/pygam/utils.py b/pygam/utils.py
index bc65ad5c..279989bc 100644
--- a/pygam/utils.py
+++ b/pygam/utils.py
@@ -64,7 +64,7 @@ def cholesky(A, sparse=True, verbose=True):  # noqa: F811
 
         if sparse:
             return L.T  # upper triangular factorization
-        return L.T.A  # upper triangular factorization
+        return L.T.toarray()  # upper triangular factorization
 
     else:
         msg = (
@@ -78,7 +78,7 @@ def cholesky(A, sparse=True, verbose=True):  # noqa: F811
             warnings.warn(msg)
 
         if sp.sparse.issparse(A):
-            A = A.A
+            A = A.toarray()
 
         try:
             L = sp.linalg.cholesky(A, lower=False)
@@ -951,10 +951,10 @@ def tensor_product(a, b, reshape=True):
         raise ValueError('both arguments must have the same number of samples')
 
     if sp.sparse.issparse(a):
-        a = a.A
+        a = a.toarray()
 
     if sp.sparse.issparse(b):
-        b = b.A
+        b = b.toarray()
 
     tensor = a[..., :, None] * b[..., None, :]