diff --git a/masque/shapes/path.py b/masque/shapes/path.py
index ff1c6e9..feb1700 100644
--- a/masque/shapes/path.py
+++ b/masque/shapes/path.py
@@ -344,8 +344,18 @@ class Path(Shape):
         bs = v[1:-1] - v[:-2] + perp[1:] - perp[:-1]
         ds = v[1:-1] - v[:-2] - perp[1:] + perp[:-1]
 
-        rp = numpy.linalg.lstsq(As, bs[:, :, None], rcond=1e-12)[0][:, 0, 0]
-        rn = numpy.linalg.lstsq(As, ds[:, :, None], rcond=1e-12)[0][:, 0, 0]
+        try:
+            # Vectorized solve for all intersections
+            # solve supports broadcasting: As (N-2, 2, 2), bs (N-2, 2, 1)
+            rp = numpy.linalg.solve(As, bs[:, :, None])[:, 0, 0]
+            rn = numpy.linalg.solve(As, ds[:, :, None])[:, 0, 0]
+        except numpy.linalg.LinAlgError:
+            # Fallback to slower lstsq if some segments are parallel (singular matrix)
+            rp = numpy.zeros(As.shape[0])
+            rn = numpy.zeros(As.shape[0])
+            for ii in range(As.shape[0]):
+                rp[ii] = numpy.linalg.lstsq(As[ii], bs[ii, :, None], rcond=1e-12)[0][0, 0]
+                rn[ii] = numpy.linalg.lstsq(As[ii], ds[ii, :, None], rcond=1e-12)[0][0, 0]
 
         intersection_p = v[:-2] + rp[:, None] * dv[:-1] + perp[:-1]
         intersection_n = v[:-2] + rn[:, None] * dv[:-1] - perp[:-1]