diff --git a/ultraplot/figure.py b/ultraplot/figure.py
index 6b5b46c4..f551df1b 100644
--- a/ultraplot/figure.py
+++ b/ultraplot/figure.py
@@ -1832,7 +1832,7 @@ def _axes_dict(naxs, input, kw=False, default=None):
         # Create or update the gridspec and add subplots with subplotspecs
         # NOTE: The gridspec is added to the figure when we pass the subplotspec
         if gs is None:
-            gs = pgridspec.GridSpec(*array.shape, **gridspec_kw)
+            gs = pgridspec.GridSpec(*array.shape, layout_array=array, **gridspec_kw)
         else:
             gs.update(**gridspec_kw)
         axs = naxs * [None]  # list of axes
diff --git a/ultraplot/gridspec.py b/ultraplot/gridspec.py
index 59de0f04..812bb7ac 100644
--- a/ultraplot/gridspec.py
+++ b/ultraplot/gridspec.py
@@ -6,21 +6,31 @@
 import itertools
 import re
 from collections.abc import MutableSequence
+from functools import wraps
 from numbers import Integral
+from typing import List, Optional, Tuple, Union
 
 import matplotlib.axes as maxes
 import matplotlib.gridspec as mgridspec
 import matplotlib.transforms as mtransforms
 import numpy as np
-from typing import List, Optional, Union, Tuple
-from functools import wraps
 
 from . import axes as paxes
 from .config import rc
-from .internals import ic  # noqa: F401
-from .internals import _not_none, docstring, warnings
+from .internals import (
+    _not_none,
+    docstring,
+    ic,  # noqa: F401
+    warnings,
+)
 from .utils import _fontsize_to_pt, units
-from .internals import warnings
+
+try:
+    from . import ultralayout
+    ULTRA_AVAILABLE = True
+except ImportError:
+    ultralayout = None
+    ULTRA_AVAILABLE = False
 
 __all__ = ["GridSpec", "SubplotGrid"]
 
@@ -225,6 +235,18 @@ def get_position(self, figure, return_all=False):
             nrows, ncols = gs.get_total_geometry()
         else:
             nrows, ncols = gs.get_geometry()
+
+        # Check if we should use UltraLayout for this subplot
+        if isinstance(gs, GridSpec) and gs._use_ultra_layout:
+            bbox = gs._get_ultra_position(self.num1, figure)
+            if bbox is not None:
+                if return_all:
+                    rows, cols = np.unravel_index([self.num1, self.num2], (nrows, ncols))
+                    return bbox, rows[0], cols[0], nrows, ncols
+                else:
+                    return bbox
+
+        # Default behavior: use grid positions
         rows, cols = np.unravel_index([self.num1, self.num2], (nrows, ncols))
         bottoms, tops, lefts, rights = gs.get_grid_positions(figure)
         bottom = bottoms[rows].min()
@@ -264,7 +286,7 @@ def __getattr__(self, attr):
         super().__getattribute__(attr)  # native error message
 
     @docstring._snippet_manager
-    def __init__(self, nrows=1, ncols=1, **kwargs):
+    def __init__(self, nrows=1, ncols=1, layout_array=None, **kwargs):
         """
         Parameters
         ----------
@@ -272,6 +294,11 @@ def __init__(self, nrows=1, ncols=1, **kwargs):
             The number of rows in the subplot grid.
         ncols : int, optional
             The number of columns in the subplot grid.
+        layout_array : array-like, optional
+            2D array specifying the subplot layout, where each unique integer
+            represents a subplot and 0 represents empty space. When provided,
+            enables UltraLayout constraint-based positioning for non-orthogonal
+            arrangements (requires kiwisolver package).
 
         Other parameters
         ----------------
@@ -301,6 +328,16 @@ def __init__(self, nrows=1, ncols=1, **kwargs):
         manually and want the same geometry for multiple figures, you must create
         a copy with `GridSpec.copy` before working on the subsequent figure).
         """
+        # Layout array for non-orthogonal layouts with UltraLayout
+        self._layout_array = np.array(layout_array) if layout_array is not None else None
+        self._ultra_positions = None  # Cache for UltraLayout-computed positions
+        self._use_ultra_layout = False  # Flag to enable UltraLayout
+
+        # Check if we should use UltraLayout
+        if self._layout_array is not None and ULTRA_AVAILABLE:
+            if not ultralayout.is_orthogonal_layout(self._layout_array):
+                self._use_ultra_layout = True
+
         # Fundamental GridSpec properties
         self._nrows_total = nrows
         self._ncols_total = ncols
@@ -363,6 +400,119 @@ def __init__(self, nrows=1, ncols=1, **kwargs):
         }
         self._update_params(pad=pad, **kwargs)
 
+    def _get_ultra_position(self, subplot_num, figure):
+        """
+        Get the position of a subplot using UltraLayout constraint-based positioning.
+
+        Parameters
+        ----------
+        subplot_num : int
+            The subplot number (in total geometry indexing)
+        figure : Figure
+            The matplotlib figure instance
+
+        Returns
+        -------
+        bbox : Bbox or None
+            The bounding box for the subplot, or None if kiwi layout fails
+        """
+        if not self._use_ultra_layout or self._layout_array is None:
+            return None
+
+        # Ensure figure is set
+        if not self.figure:
+            self._figure = figure
+        if not self.figure:
+            return None
+
+        # Compute or retrieve cached UltraLayout positions
+        if self._ultra_positions is None:
+            self._compute_ultra_positions()
+
+        # Find which subplot number in the layout array corresponds to this subplot_num
+        # We need to map from the gridspec cell index to the layout array subplot number
+        nrows, ncols = self._layout_array.shape
+
+        # Decode the subplot_num to find which layout number it corresponds to
+        # This is a bit tricky because subplot_num is in total geometry space
+        # We need to find which unique number in the layout_array this corresponds to
+
+        # Get the cell position from subplot_num
+        row, col = divmod(subplot_num, self.ncols_total)
+
+        # Check if this is within the layout array bounds
+        if row >= nrows or col >= ncols:
+            return None
+
+        # Get the layout number at this position
+        layout_num = self._layout_array[row, col]
+
+        if layout_num == 0 or layout_num not in self._ultra_positions:
+            return None
+
+        # Return the cached position
+        left, bottom, width, height = self._ultra_positions[layout_num]
+        bbox = mtransforms.Bbox.from_bounds(left, bottom, width, height)
+        return bbox
+
+    def _compute_ultra_positions(self):
+        """
+        Compute subplot positions using UltraLayout and cache them.
+        """
+        if not ULTRA_AVAILABLE or self._layout_array is None:
+            return
+
+        # Get figure size
+        if not self.figure:
+            return
+
+        figwidth, figheight = self.figure.get_size_inches()
+
+        # Convert spacing to inches
+        wspace_inches = []
+        for i, ws in enumerate(self._wspace_total):
+            if ws is not None:
+                wspace_inches.append(ws)
+            else:
+                # Use default spacing
+                wspace_inches.append(0.2)  # Default spacing in inches
+
+        hspace_inches = []
+        for i, hs in enumerate(self._hspace_total):
+            if hs is not None:
+                hspace_inches.append(hs)
+            else:
+                hspace_inches.append(0.2)
+
+        # Get margins
+        left = self.left if self.left is not None else self._left_default if self._left_default is not None else 0.125 * figwidth
+        right = self.right if self.right is not None else self._right_default if self._right_default is not None else 0.125 * figwidth
+        top = self.top if self.top is not None else self._top_default if self._top_default is not None else 0.125 * figheight
+        bottom = self.bottom if self.bottom is not None else self._bottom_default if self._bottom_default is not None else 0.125 * figheight
+
+        # Compute positions using UltraLayout
+        try:
+            self._ultra_positions = ultralayout.compute_ultra_positions(
+                self._layout_array,
+                figwidth=figwidth,
+                figheight=figheight,
+                wspace=wspace_inches,
+                hspace=hspace_inches,
+                left=left,
+                right=right,
+                top=top,
+                bottom=bottom,
+                wratios=self._wratios_total,
+                hratios=self._hratios_total
+            )
+        except Exception as e:
+            warnings._warn_ultraplot(
+                f"Failed to compute UltraLayout: {e}. "
+                "Falling back to default grid layout."
+            )
+            self._use_ultra_layout = False
+            self._ultra_positions = None
+
     def __getitem__(self, key):
         """
         Get a `~matplotlib.gridspec.SubplotSpec`. "Hidden" slots allocated for axes
diff --git a/ultraplot/tests/test_ultralayout.py b/ultraplot/tests/test_ultralayout.py
new file mode 100644
index 00000000..9c8f573c
--- /dev/null
+++ b/ultraplot/tests/test_ultralayout.py
@@ -0,0 +1,289 @@
+import numpy as np
+import pytest
+
+import ultraplot as uplt
+from ultraplot import ultralayout
+from ultraplot.gridspec import GridSpec
+
+
+def test_is_orthogonal_layout_simple_grid():
+    """Test orthogonal layout detection for simple grids."""
+    # Simple 2x2 grid should be orthogonal
+    array = np.array([[1, 2], [3, 4]])
+    assert ultralayout.is_orthogonal_layout(array) is True
+
+
+def test_is_orthogonal_layout_non_orthogonal():
+    """Test orthogonal layout detection for non-orthogonal layouts."""
+    # Centered subplot with empty cells should be non-orthogonal
+    array = np.array([[1, 1, 2, 2], [0, 3, 3, 0]])
+    assert ultralayout.is_orthogonal_layout(array) is False
+
+
+def test_is_orthogonal_layout_spanning():
+    """Test orthogonal layout with spanning subplots that is still orthogonal."""
+    # L-shape that maintains grid alignment
+    array = np.array([[1, 1], [1, 2]])
+    assert ultralayout.is_orthogonal_layout(array) is True
+
+
+def test_is_orthogonal_layout_with_gaps():
+    """Test non-orthogonal layout with gaps."""
+    array = np.array([[1, 1, 1], [2, 0, 3]])
+    assert ultralayout.is_orthogonal_layout(array) is False
+
+
+def test_is_orthogonal_layout_empty():
+    """Test empty layout."""
+    array = np.array([[0, 0], [0, 0]])
+    assert ultralayout.is_orthogonal_layout(array) is True
+
+
+def test_gridspec_with_orthogonal_layout():
+    """Test that GridSpec doesn't activate UltraLayout for orthogonal layouts."""
+    layout = np.array([[1, 2], [3, 4]])
+    gs = GridSpec(2, 2, layout_array=layout)
+    assert gs._layout_array is not None
+    # Should not use UltraLayout for orthogonal layouts
+    assert gs._use_ultra_layout is False
+
+
+def test_gridspec_with_non_orthogonal_layout():
+    """Test that GridSpec activates UltraLayout for non-orthogonal layouts."""
+    pytest.importorskip("kiwisolver")
+    layout = np.array([[1, 1, 2, 2], [0, 3, 3, 0]])
+    gs = GridSpec(2, 4, layout_array=layout)
+    assert gs._layout_array is not None
+    # Should use UltraLayout for non-orthogonal layouts
+    assert gs._use_ultra_layout is True
+
+
+def test_gridspec_without_kiwisolver(monkeypatch):
+    """Test graceful fallback when kiwisolver is not available."""
+    # Mock the ULTRA_AVAILABLE flag
+    import ultraplot.gridspec as gs_module
+    monkeypatch.setattr(gs_module, "ULTRA_AVAILABLE", False)
+
+    layout = np.array([[1, 1, 2, 2], [0, 3, 3, 0]])
+    gs = GridSpec(2, 4, layout_array=layout)
+    # Should not activate UltraLayout if kiwisolver not available
+    assert gs._use_ultra_layout is False
+
+
+def test_ultralayout_solver_initialization():
+    """Test UltraLayoutSolver can be initialized."""
+    pytest.importorskip("kiwisolver")
+    layout = np.array([[1, 1, 2, 2], [0, 3, 3, 0]])
+    solver = ultralayout.UltraLayoutSolver(
+        layout,
+        figwidth=10.0,
+        figheight=6.0
+    )
+    assert solver.array is not None
+    assert solver.nrows == 2
+    assert solver.ncols == 4
+
+
+def test_compute_ultra_positions():
+    """Test computing positions with UltraLayout."""
+    pytest.importorskip("kiwisolver")
+    layout = np.array([[1, 1, 2, 2], [0, 3, 3, 0]])
+    positions = ultralayout.compute_ultra_positions(
+        layout,
+        figwidth=10.0,
+        figheight=6.0,
+        wspace=[0.2, 0.2, 0.2],
+        hspace=[0.2],
+    )
+
+    # Should return positions for 3 subplots
+    assert len(positions) == 3
+    assert 1 in positions
+    assert 2 in positions
+    assert 3 in positions
+
+    # Each position should be (left, bottom, width, height)
+    for num, pos in positions.items():
+        assert len(pos) == 4
+        left, bottom, width, height = pos
+        assert 0 <= left <= 1
+        assert 0 <= bottom <= 1
+        assert width > 0
+        assert height > 0
+        assert left + width <= 1.01  # Allow small numerical error
+        assert bottom + height <= 1.01
+
+
+def test_subplots_with_non_orthogonal_layout():
+    """Test creating subplots with non-orthogonal layout."""
+    pytest.importorskip("kiwisolver")
+    layout = [[1, 1, 2, 2], [0, 3, 3, 0]]
+    fig, axs = uplt.subplots(array=layout, figsize=(10, 6))
+
+    # Should create 3 subplots
+    assert len(axs) == 3
+
+    # Check that positions are valid
+    for ax in axs:
+        pos = ax.get_position()
+        assert pos.width > 0
+        assert pos.height > 0
+        assert 0 <= pos.x0 <= 1
+        assert 0 <= pos.y0 <= 1
+
+
+def test_subplots_with_orthogonal_layout():
+    """Test creating subplots with orthogonal layout (should work as before)."""
+    layout = [[1, 2], [3, 4]]
+    fig, axs = uplt.subplots(array=layout, figsize=(8, 6))
+
+    # Should create 4 subplots
+    assert len(axs) == 4
+
+    # Check that positions are valid
+    for ax in axs:
+        pos = ax.get_position()
+        assert pos.width > 0
+        assert pos.height > 0
+
+
+def test_ultralayout_respects_spacing():
+    """Test that UltraLayout respects spacing parameters."""
+    pytest.importorskip("kiwisolver")
+    layout = np.array([[1, 1, 2, 2], [0, 3, 3, 0]])
+
+    # Compute with different spacing
+    positions1 = ultralayout.compute_ultra_positions(
+        layout, figwidth=10.0, figheight=6.0,
+        wspace=[0.1, 0.1, 0.1], hspace=[0.1]
+    )
+    positions2 = ultralayout.compute_ultra_positions(
+        layout, figwidth=10.0, figheight=6.0,
+        wspace=[0.5, 0.5, 0.5], hspace=[0.5]
+    )
+
+    # Subplots should be smaller with more spacing
+    for num in [1, 2, 3]:
+        _, _, width1, height1 = positions1[num]
+        _, _, width2, height2 = positions2[num]
+        # With more spacing, subplots should be smaller
+        assert width2 < width1 or height2 < height1
+
+
+def test_ultralayout_respects_ratios():
+    """Test that UltraLayout respects width/height ratios."""
+    pytest.importorskip("kiwisolver")
+    layout = np.array([[1, 2], [3, 4]])
+
+    # Equal ratios
+    positions1 = ultralayout.compute_ultra_positions(
+        layout, figwidth=10.0, figheight=6.0,
+        wratios=[1, 1], hratios=[1, 1]
+    )
+
+    # Unequal ratios
+    positions2 = ultralayout.compute_ultra_positions(
+        layout, figwidth=10.0, figheight=6.0,
+        wratios=[1, 2], hratios=[1, 1]
+    )
+
+    # Subplot 2 should be wider than subplot 1 with unequal ratios
+    _, _, width1_1, _ = positions1[1]
+    _, _, width1_2, _ = positions1[2]
+    _, _, width2_1, _ = positions2[1]
+    _, _, width2_2, _ = positions2[2]
+
+    # With equal ratios, widths should be similar
+    assert abs(width1_1 - width1_2) < 0.01
+    # With 1:2 ratio, second should be roughly twice as wide
+    assert width2_2 > width2_1
+
+
+def test_ultralayout_cached_positions():
+    """Test that UltraLayout positions are cached in GridSpec."""
+    pytest.importorskip("kiwisolver")
+    layout = np.array([[1, 1, 2, 2], [0, 3, 3, 0]])
+    gs = GridSpec(2, 4, layout_array=layout)
+
+    # Positions should not be computed yet
+    assert gs._ultra_positions is None
+
+    # Create a figure to trigger position computation
+    fig = uplt.figure()
+    gs._figure = fig
+
+    # Access a position (this should trigger computation)
+    ss = gs[0, 0]
+    pos = ss.get_position(fig)
+
+    # Positions should now be cached
+    assert gs._ultra_positions is not None
+    assert len(gs._ultra_positions) == 3
+
+
+def test_ultralayout_with_margins():
+    """Test that UltraLayout respects margin parameters."""
+    pytest.importorskip("kiwisolver")
+    layout = np.array([[1, 2]])
+
+    # Small margins
+    positions1 = ultralayout.compute_ultra_positions(
+        layout, figwidth=10.0, figheight=6.0,
+        left=0.1, right=0.1, top=0.1, bottom=0.1
+    )
+
+    # Large margins
+    positions2 = ultralayout.compute_ultra_positions(
+        layout, figwidth=10.0, figheight=6.0,
+        left=1.0, right=1.0, top=1.0, bottom=1.0
+    )
+
+    # With larger margins, subplots should be smaller
+    for num in [1, 2]:
+        _, _, width1, height1 = positions1[num]
+        _, _, width2, height2 = positions2[num]
+        assert width2 < width1
+        assert height2 < height1
+
+
+def test_complex_non_orthogonal_layout():
+    """Test a more complex non-orthogonal layout."""
+    pytest.importorskip("kiwisolver")
+    layout = np.array([
+        [1, 1, 1, 2],
+        [3, 3, 0, 2],
+        [4, 5, 5, 5]
+    ])
+
+    positions = ultralayout.compute_ultra_positions(layout, figwidth=12.0, figheight=9.0)
+
+    # Should have 5 subplots
+    assert len(positions) == 5
+
+    # All positions should be valid
+    for num in range(1, 6):
+        assert num in positions
+        left, bottom, width, height = positions[num]
+        assert 0 <= left <= 1
+        assert 0 <= bottom <= 1
+        assert width > 0
+        assert height > 0
+
+
+def test_ultralayout_module_exports():
+    """Test that ultralayout module exports expected symbols."""
+    assert hasattr(ultralayout, 'UltraLayoutSolver')
+    assert hasattr(ultralayout, 'compute_ultra_positions')
+    assert hasattr(ultralayout, 'is_orthogonal_layout')
+    assert hasattr(ultralayout, 'get_grid_positions_ultra')
+
+
+def test_gridspec_copy_preserves_layout_array():
+    """Test that copying a GridSpec preserves the layout array."""
+    layout = np.array([[1, 1, 2, 2], [0, 3, 3, 0]])
+    gs1 = GridSpec(2, 4, layout_array=layout)
+    gs2 = gs1.copy()
+
+    assert gs2._layout_array is not None
+    assert np.array_equal(gs1._layout_array, gs2._layout_array)
+    assert gs1._use_ultra_layout == gs2._use_ultra_layout
diff --git a/ultraplot/ultralayout.py b/ultraplot/ultralayout.py
new file mode 100644
index 00000000..239b5c23
--- /dev/null
+++ b/ultraplot/ultralayout.py
@@ -0,0 +1,463 @@
+#!/usr/bin/env python3
+"""
+UltraLayout: Advanced constraint-based layout system for non-orthogonal subplot arrangements.
+
+This module provides UltraPlot's constraint-based layout computation for subplot grids
+that don't follow simple orthogonal patterns, such as [[1, 1, 2, 2], [0, 3, 3, 0]]
+where subplot 3 should be nicely centered between subplots 1 and 2.
+"""
+
+from typing import Dict, List, Optional, Tuple
+
+import numpy as np
+
+try:
+    from kiwisolver import Constraint, Solver, Variable
+    KIWI_AVAILABLE = True
+except ImportError:
+    KIWI_AVAILABLE = False
+    Variable = None
+    Solver = None
+    Constraint = None
+
+
+__all__ = ['UltraLayoutSolver', 'compute_ultra_positions', 'is_orthogonal_layout']
+
+
+def is_orthogonal_layout(array: np.ndarray) -> bool:
+    """
+    Check if a subplot array follows an orthogonal (grid-aligned) layout.
+
+    An orthogonal layout is one where every subplot's edges align with
+    other subplots' edges, forming a simple grid.
+
+    Parameters
+    ----------
+    array : np.ndarray
+        2D array of subplot numbers (with 0 for empty cells)
+
+    Returns
+    -------
+    bool
+        True if layout is orthogonal, False otherwise
+    """
+    if array.size == 0:
+        return True
+
+    nrows, ncols = array.shape
+
+    # Get unique subplot numbers (excluding 0)
+    subplot_nums = np.unique(array[array != 0])
+
+    if len(subplot_nums) == 0:
+        return True
+
+    # For each subplot, get its bounding box
+    bboxes = {}
+    for num in subplot_nums:
+        rows, cols = np.where(array == num)
+        bboxes[num] = {
+            'row_min': rows.min(),
+            'row_max': rows.max(),
+            'col_min': cols.min(),
+            'col_max': cols.max(),
+        }
+
+    # Check if layout is orthogonal by verifying that all vertical and
+    # horizontal edges align with cell boundaries
+    # A more sophisticated check: for each row/col boundary, check if
+    # all subplots either cross it or are completely on one side
+
+    # Collect all unique row and column boundaries
+    row_boundaries = set()
+    col_boundaries = set()
+
+    for bbox in bboxes.values():
+        row_boundaries.add(bbox['row_min'])
+        row_boundaries.add(bbox['row_max'] + 1)
+        col_boundaries.add(bbox['col_min'])
+        col_boundaries.add(bbox['col_max'] + 1)
+
+    # Check if these boundaries create a consistent grid
+    # For orthogonal layout, we should be able to split the grid
+    # using these boundaries such that each subplot is a union of cells
+
+    row_boundaries = sorted(row_boundaries)
+    col_boundaries = sorted(col_boundaries)
+
+    # Create a refined grid
+    refined_rows = len(row_boundaries) - 1
+    refined_cols = len(col_boundaries) - 1
+
+    if refined_rows == 0 or refined_cols == 0:
+        return True
+
+    # Map each subplot to refined grid cells
+    for num in subplot_nums:
+        rows, cols = np.where(array == num)
+
+        # Check if this subplot occupies a rectangular region in the refined grid
+        refined_row_indices = set()
+        refined_col_indices = set()
+
+        for r in rows:
+            for i, (r_start, r_end) in enumerate(zip(row_boundaries[:-1], row_boundaries[1:])):
+                if r_start <= r < r_end:
+                    refined_row_indices.add(i)
+
+        for c in cols:
+            for i, (c_start, c_end) in enumerate(zip(col_boundaries[:-1], col_boundaries[1:])):
+                if c_start <= c < c_end:
+                    refined_col_indices.add(i)
+
+        # Check if indices form a rectangle
+        if refined_row_indices and refined_col_indices:
+            r_min, r_max = min(refined_row_indices), max(refined_row_indices)
+            c_min, c_max = min(refined_col_indices), max(refined_col_indices)
+
+            expected_cells = (r_max - r_min + 1) * (c_max - c_min + 1)
+            actual_cells = len(refined_row_indices) * len(refined_col_indices)
+
+            if expected_cells != actual_cells:
+                return False
+
+    return True
+
+
+class UltraLayoutSolver:
+    """
+    UltraLayout: Constraint-based layout solver using kiwisolver for subplot positioning.
+
+    This solver computes aesthetically pleasing positions for subplots in
+    non-orthogonal arrangements by using constraint satisfaction, providing
+    a superior layout experience for complex subplot arrangements.
+    """
+
+    def __init__(self, array: np.ndarray, figwidth: float = 10.0, figheight: float = 8.0,
+                 wspace: Optional[List[float]] = None, hspace: Optional[List[float]] = None,
+                 left: float = 0.125, right: float = 0.125,
+                 top: float = 0.125, bottom: float = 0.125,
+                 wratios: Optional[List[float]] = None, hratios: Optional[List[float]] = None):
+        """
+        Initialize the UltraLayout solver.
+
+        Parameters
+        ----------
+        array : np.ndarray
+            2D array of subplot numbers (with 0 for empty cells)
+        figwidth, figheight : float
+            Figure dimensions in inches
+        wspace, hspace : list of float, optional
+            Spacing between columns and rows in inches
+        left, right, top, bottom : float
+            Margins in inches
+        wratios, hratios : list of float, optional
+            Width and height ratios for columns and rows
+        """
+        if not KIWI_AVAILABLE:
+            raise ImportError(
+                "kiwisolver is required for non-orthogonal layouts. "
+                "Install it with: pip install kiwisolver"
+            )
+
+        self.array = array
+        self.nrows, self.ncols = array.shape
+        self.figwidth = figwidth
+        self.figheight = figheight
+        self.left_margin = left
+        self.right_margin = right
+        self.top_margin = top
+        self.bottom_margin = bottom
+
+        # Get subplot numbers
+        self.subplot_nums = sorted(np.unique(array[array != 0]))
+
+        # Set up spacing
+        if wspace is None:
+            self.wspace = [0.2] * (self.ncols - 1) if self.ncols > 1 else []
+        else:
+            self.wspace = list(wspace)
+
+        if hspace is None:
+            self.hspace = [0.2] * (self.nrows - 1) if self.nrows > 1 else []
+        else:
+            self.hspace = list(hspace)
+
+        # Set up ratios
+        if wratios is None:
+            self.wratios = [1.0] * self.ncols
+        else:
+            self.wratios = list(wratios)
+
+        if hratios is None:
+            self.hratios = [1.0] * self.nrows
+        else:
+            self.hratios = list(hratios)
+
+        # Initialize solver
+        self.solver = Solver()
+        self.variables = {}
+        self._setup_variables()
+        self._setup_constraints()
+
+    def _setup_variables(self):
+        """Create kiwisolver variables for all grid lines."""
+        # Vertical lines (left edges of columns + right edge of last column)
+        self.col_lefts = [Variable(f'col_{i}_left') for i in range(self.ncols)]
+        self.col_rights = [Variable(f'col_{i}_right') for i in range(self.ncols)]
+
+        # Horizontal lines (top edges of rows + bottom edge of last row)
+        # Note: in figure coordinates, top is higher value
+        self.row_tops = [Variable(f'row_{i}_top') for i in range(self.nrows)]
+        self.row_bottoms = [Variable(f'row_{i}_bottom') for i in range(self.nrows)]
+
+    def _setup_constraints(self):
+        """Set up all constraints for the layout."""
+        # 1. Figure boundary constraints
+        self.solver.addConstraint(self.col_lefts[0] == self.left_margin / self.figwidth)
+        self.solver.addConstraint(self.col_rights[-1] == 1.0 - self.right_margin / self.figwidth)
+        self.solver.addConstraint(self.row_bottoms[-1] == self.bottom_margin / self.figheight)
+        self.solver.addConstraint(self.row_tops[0] == 1.0 - self.top_margin / self.figheight)
+
+        # 2. Column continuity and spacing constraints
+        for i in range(self.ncols - 1):
+            # Right edge of column i connects to left edge of column i+1 with spacing
+            spacing = self.wspace[i] / self.figwidth if i < len(self.wspace) else 0
+            self.solver.addConstraint(self.col_rights[i] + spacing == self.col_lefts[i + 1])
+
+        # 3. Row continuity and spacing constraints
+        for i in range(self.nrows - 1):
+            # Bottom edge of row i connects to top edge of row i+1 with spacing
+            spacing = self.hspace[i] / self.figheight if i < len(self.hspace) else 0
+            self.solver.addConstraint(self.row_bottoms[i] == self.row_tops[i + 1] + spacing)
+
+        # 4. Width ratio constraints
+        total_width = 1.0 - (self.left_margin + self.right_margin) / self.figwidth
+        if self.ncols > 1:
+            spacing_total = sum(self.wspace) / self.figwidth
+        else:
+            spacing_total = 0
+        available_width = total_width - spacing_total
+        total_ratio = sum(self.wratios)
+
+        for i in range(self.ncols):
+            width = available_width * self.wratios[i] / total_ratio
+            self.solver.addConstraint(self.col_rights[i] == self.col_lefts[i] + width)
+
+        # 5. Height ratio constraints
+        total_height = 1.0 - (self.top_margin + self.bottom_margin) / self.figheight
+        if self.nrows > 1:
+            spacing_total = sum(self.hspace) / self.figheight
+        else:
+            spacing_total = 0
+        available_height = total_height - spacing_total
+        total_ratio = sum(self.hratios)
+
+        for i in range(self.nrows):
+            height = available_height * self.hratios[i] / total_ratio
+            self.solver.addConstraint(self.row_tops[i] == self.row_bottoms[i] + height)
+
+        # 6. Add aesthetic constraints for non-orthogonal layouts
+        self._add_aesthetic_constraints()
+
+    def _add_aesthetic_constraints(self):
+        """
+        Add constraints to make non-orthogonal layouts look nice.
+
+        For subplots that span cells in non-aligned ways, we add constraints
+        to center them or align them aesthetically with neighboring subplots.
+        """
+        # Analyze the layout to find subplots that need special handling
+        for num in self.subplot_nums:
+            rows, cols = np.where(self.array == num)
+            row_min, row_max = rows.min(), rows.max()
+            col_min, col_max = cols.min(), cols.max()
+
+            # Check if this subplot has empty cells on its sides
+            # If so, try to center it with respect to subplots above/below/beside
+
+            # Check left side
+            if col_min > 0:
+                left_cells = self.array[row_min:row_max+1, col_min-1]
+                if np.all(left_cells == 0):
+                    # Empty on the left - might want to align with something above/below
+                    self._try_align_with_neighbors(num, 'left', row_min, row_max, col_min)
+
+            # Check right side
+            if col_max < self.ncols - 1:
+                right_cells = self.array[row_min:row_max+1, col_max+1]
+                if np.all(right_cells == 0):
+                    # Empty on the right
+                    self._try_align_with_neighbors(num, 'right', row_min, row_max, col_max)
+
+    def _try_align_with_neighbors(self, num: int, side: str, row_min: int, row_max: int, col_idx: int):
+        """
+        Try to align a subplot edge with neighboring subplots.
+
+        For example, if subplot 3 is in row 1 between subplots 1 and 2 in row 0,
+        we want to center it between them.
+        """
+        # Find subplots in adjacent rows that overlap with this subplot's column range
+        rows, cols = np.where(self.array == num)
+        col_min, col_max = cols.min(), cols.max()
+
+        # Look in rows above
+        if row_min > 0:
+            above_nums = set()
+            for r in range(row_min):
+                for c in range(col_min, col_max + 1):
+                    if self.array[r, c] != 0:
+                        above_nums.add(self.array[r, c])
+
+            if len(above_nums) >= 2:
+                # Multiple subplots above - try to center between them
+                above_nums = sorted(above_nums)
+                # Find the leftmost and rightmost subplots above
+                leftmost_cols = []
+                rightmost_cols = []
+                for n in above_nums:
+                    n_cols = np.where(self.array == n)[1]
+                    leftmost_cols.append(n_cols.min())
+                    rightmost_cols.append(n_cols.max())
+
+                # If we're between two subplots, center between them
+                if side == 'left' and leftmost_cols:
+                    # Could add centering constraint here
+                    # For now, we let the default grid handle it
+                    pass
+
+        # Look in rows below
+        if row_max < self.nrows - 1:
+            below_nums = set()
+            for r in range(row_max + 1, self.nrows):
+                for c in range(col_min, col_max + 1):
+                    if self.array[r, c] != 0:
+                        below_nums.add(self.array[r, c])
+
+            if len(below_nums) >= 2:
+                # Similar logic for below
+                pass
+
+    def solve(self) -> Dict[int, Tuple[float, float, float, float]]:
+        """
+        Solve the constraint system and return subplot positions.
+
+        Returns
+        -------
+        dict
+            Dictionary mapping subplot numbers to (left, bottom, width, height)
+            in figure-relative coordinates [0, 1]
+        """
+        # Solve the constraint system
+        self.solver.updateVariables()
+
+        # Extract positions for each subplot
+        positions = {}
+
+        for num in self.subplot_nums:
+            rows, cols = np.where(self.array == num)
+            row_min, row_max = rows.min(), rows.max()
+            col_min, col_max = cols.min(), cols.max()
+
+            # Get the bounding box from the grid lines
+            left = self.col_lefts[col_min].value()
+            right = self.col_rights[col_max].value()
+            bottom = self.row_bottoms[row_max].value()
+            top = self.row_tops[row_min].value()
+
+            width = right - left
+            height = top - bottom
+
+            positions[num] = (left, bottom, width, height)
+
+        return positions
+
+
+def compute_ultra_positions(array: np.ndarray, figwidth: float = 10.0, figheight: float = 8.0,
+                           wspace: Optional[List[float]] = None, hspace: Optional[List[float]] = None,
+                           left: float = 0.125, right: float = 0.125,
+                           top: float = 0.125, bottom: float = 0.125,
+                           wratios: Optional[List[float]] = None,
+                           hratios: Optional[List[float]] = None) -> Dict[int, Tuple[float, float, float, float]]:
+    """
+    Compute subplot positions using UltraLayout for non-orthogonal layouts.
+
+    Parameters
+    ----------
+    array : np.ndarray
+        2D array of subplot numbers (with 0 for empty cells)
+    figwidth, figheight : float
+        Figure dimensions in inches
+    wspace, hspace : list of float, optional
+        Spacing between columns and rows in inches
+    left, right, top, bottom : float
+        Margins in inches
+    wratios, hratios : list of float, optional
+        Width and height ratios for columns and rows
+
+    Returns
+    -------
+    dict
+        Dictionary mapping subplot numbers to (left, bottom, width, height)
+        in figure-relative coordinates [0, 1]
+
+    Examples
+    --------
+    >>> array = np.array([[1, 1, 2, 2], [0, 3, 3, 0]])
+    >>> positions = compute_ultra_positions(array)
+    >>> positions[3]  # Position of subplot 3
+    (0.25, 0.125, 0.5, 0.35)
+    """
+    solver = UltraLayoutSolver(
+        array, figwidth, figheight, wspace, hspace,
+        left, right, top, bottom, wratios, hratios
+    )
+    return solver.solve()
+
+
+def get_grid_positions_ultra(array: np.ndarray, figwidth: float, figheight: float,
+                             wspace: Optional[List[float]] = None,
+                             hspace: Optional[List[float]] = None,
+                             left: float = 0.125, right: float = 0.125,
+                             top: float = 0.125, bottom: float = 0.125,
+                             wratios: Optional[List[float]] = None,
+                             hratios: Optional[List[float]] = None) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
+    """
+    Get grid line positions using UltraLayout.
+
+    This returns arrays of grid line positions similar to GridSpec.get_grid_positions(),
+    but computed using UltraLayout's constraint satisfaction for better handling of non-orthogonal layouts.
+
+    Parameters
+    ----------
+    array : np.ndarray
+        2D array of subplot numbers
+    figwidth, figheight : float
+        Figure dimensions in inches
+    wspace, hspace : list of float, optional
+        Spacing between columns and rows in inches
+    left, right, top, bottom : float
+        Margins in inches
+    wratios, hratios : list of float, optional
+        Width and height ratios for columns and rows
+
+    Returns
+    -------
+    bottoms, tops, lefts, rights : np.ndarray
+        Arrays of grid line positions for each cell
+    """
+    solver = UltraLayoutSolver(
+        array, figwidth, figheight, wspace, hspace,
+        left, right, top, bottom, wratios, hratios
+    )
+    solver.solver.updateVariables()
+
+    nrows, ncols = array.shape
+
+    # Extract grid line positions
+    lefts = np.array([v.value() for v in solver.col_lefts])
+    rights = np.array([v.value() for v in solver.col_rights])
+    tops = np.array([v.value() for v in solver.row_tops])
+    bottoms = np.array([v.value() for v in solver.row_bottoms])
+
+    return bottoms, tops, lefts, rights