Return to Answer

SpatiaLite implements a new way to query the spatial index, through the SpatialIndex table. It is a virtual table wrapping the rectangles tree, and it is presented in: https://gaia-gis.it/fossil/libspatialite/wiki?name=SpatialIndex.
But neither the GeoPackage standard (nor its GDAL implementation) includes it.

By querying the spatial index, it is very quick to preselect features whose geometries are candidates for the intersection to return true.

SpatiaLite implements a new way to query the spatial index, through the SpatialIndex table. It is a virtual table wrapping the rectangles tree, and it is presented in: https://gaia-gis.it/fossil/libspatialite/wiki?name=SpatialIndex.
But neither the GeoPackage standard (nor its GDAL implementation) includes it. Instead, the spatial index is wrapped in SQLite's own rtree virtual tables, which are documented in: https://www.sqlite.org/rtree.html.

import geopandas as gpd
import shapely.geometry
import numpy as np
import time

np.random.seed(42)

elapsed_with_rtree = []
elapsed_without_rtree = []

ns = [8,16, 24, 32, 40, 48, 56, 64]

for n in ns:

    p = gpd.GeoSeries(
        [shapely.geometry.box(j, i, j + 1, i + 1) for i in range(n) for j in range(n)]
    )

    q = gpd.GeoSeries(
        [
            shapely.geometry.Point(e)
            for e in np.random.uniform(low=0, high=n, size=[n * n * 10, 2])
        ]
    ).buffer(distance=0.1)
    
    !rm -f united.gpkg
    p.to_file("united.gpkg",layer="p")
    q.to_file("united.gpkg",layer="q")
    
    query = """
        SELECT ST_Difference(
          p.geom,
          (SELECT ST_UNION(geom)
           FROM q
           WHERE ST_Intersects(p.geom,geom)
           )
        ) AS geom FROM p;
    """
    query = ' '.join(query.split())
    
    rtree_query = """
        SELECT ST_Difference(
          p.geom,
          (SELECT ST_UNION(geom)
           FROM q
           WHERE ROWID IN(
             SELECT id
             FROM rtree_q_geom
             WHERE minx <= MbrMaxX(p.geom)
               AND maxx >= MbrMinX(p.geom)
               AND miny <= MbrMaxY(p.geom)
               AND maxy >= MbrMinY(p.geom)
             )
           AND ST_Intersects(p.geom, geom)
           )
        ) AS geom FROM p;
    """
    rtree_query = ' '.join(rtree_query.split())
    
    t0 = time.time()
    subprocess.run(['ogrinfo', '-q', '-sql', query, 'united.gpkg'], stdout=subprocess.DEVNULL)
    t1 = time.time()
    subprocess.run(['ogrinfo', '-q', '-sql', rtree_query, 'united.gpkg'], stdout=subprocess.DEVNULL)
    t2 = time.time()

    elapsed_without_rtree.append(t1-t0)
    elapsed_with_rtree.append(t2-t1)


import matplotlib.pyplot as plt

plt.figure(figsize=(8,4))

plt.scatter([e**2 for e in ns], elapsed_without_rtree, c='k')
plt.plot([e**2 for e in ns], elapsed_without_rtree, c='r', label="without rtree")

plt.scatter([e**2 for e in ns], elapsed_with_rtree, c='k')
plt.plot([e**2 for e in ns], elapsed_with_rtree, c='b', label="with rtree")

plt.xlabel("Number of geometries in p")
plt.ylabel("ogrinfo execution time")

plt.legend()

import geopandas as gpd
import shapely.geometry
import numpy as np
import time

np.random.seed(42)

elapsed_with_rtree = []
elapsed_without_rtree = []

ns = [8, 16, 24, 32, 40, 48, 56, 64]

for n in ns:

    p = gpd.GeoSeries(
        [shapely.geometry.box(j, i, j + 1, i + 1) for i in range(n) for j in range(n)]
    )

    q = gpd.GeoSeries(
        [
            shapely.geometry.Point(e)
            for e in np.random.uniform(low=0, high=n, size=[n * n * 10, 2])
        ]
    ).buffer(distance=0.1)
    
    !rm -f united.gpkg
    p.to_file("united.gpkg", layer="p")
    q.to_file("united.gpkg", layer="q")
    
    query = """
        SELECT ST_Difference(
          p.geom,
          (SELECT ST_UNION(geom)
           FROM q
           WHERE ST_Intersects(p.geom,geom)
           )
        ) AS geom FROM p;
    """
    query = ' '.join(query.split())
    
    rtree_query = """
        SELECT ST_Difference(
          p.geom,
          (SELECT ST_UNION(geom)
           FROM q
           WHERE ROWID IN(
             SELECT id
             FROM rtree_q_geom
             WHERE minx <= MbrMaxX(p.geom)
               AND maxx >= MbrMinX(p.geom)
               AND miny <= MbrMaxY(p.geom)
               AND maxy >= MbrMinY(p.geom)
             )
           AND ST_Intersects(p.geom, geom)
           )
        ) AS geom FROM p;
    """
    rtree_query = ' '.join(rtree_query.split())
    
    t0 = time.time()
    subprocess.run(['ogrinfo', '-q', '-sql', query, 'united.gpkg'], stdout=subprocess.DEVNULL)
    t1 = time.time()
    subprocess.run(['ogrinfo', '-q', '-sql', rtree_query, 'united.gpkg'], stdout=subprocess.DEVNULL)
    t2 = time.time()

    elapsed_without_rtree.append(t1-t0)
    elapsed_with_rtree.append(t2-t1)


import matplotlib.pyplot as plt

plt.figure(figsize=(8,4))

plt.scatter([e**2 for e in ns], elapsed_without_rtree, c='k')
plt.plot([e**2 for e in ns], elapsed_without_rtree, c='r', label="without rtree")

plt.scatter([e**2 for e in ns], elapsed_with_rtree, c='k')
plt.plot([e**2 for e in ns], elapsed_with_rtree, c='b', label="with rtree")

plt.xlabel("Number of geometries in p")
plt.ylabel("ogrinfo execution time")

plt.legend()

SELECT ST_Difference(p.geom,
 (SELECT ST_Union(geom)
  FROM q
  WHERE ST_Intersects(p.geom, geom)
  AND p.ROWID IN (
   SELECT ROWID 
    FROM SpatialIndex
   WHERE f_table_name = 'p' 
   AND search_frame = geom)
 )
) AS geom FROM p

Spatial indexes are trees where each node represent a rectangle that belongs within its parent node.

Regarding the SpatialIndex table, it is a virtual table wrapping the rectangles tree, and it is presented in: https://gaia-gis.it/fossil/libspatialite/wiki?name=SpatialIndex.

Regarding the ROWID attribute, it is documented in https://www.sqlite.org/lang_createtable.html#rowid. In the GeoPackage standard, the Feature ID attribute refers to it.

SELECT ST_Difference(
  p.geom,
  (SELECT ST_UNION(geom)
   FROM q
   WHERE ROWID IN(
     SELECT id
     FROM rtree_q_geom
     WHERE minx <= MbrMaxX(p.geom)
       AND maxx >= MbrMinX(p.geom)
       AND miny <= MbrMaxY(p.geom)
       AND maxy >= MbrMinY(p.geom)
     )
   AND ST_Intersects(p.geom, geom)
   )
) AS geom FROM p;

As Pieter answer, there are many ways to query the spatial index.

Spatial indexes are trees where each node represent a rectangle that belongs within its parent node. A brief, outdated but useful introduction was described in: http://www.gaia-gis.it/gaia-sins/spatialite-cookbook/html/rtree.html

SpatiaLite implements a new way to query the spatial index, through the SpatialIndex table. It is a virtual table wrapping the rectangles tree, and it is presented in: https://gaia-gis.it/fossil/libspatialite/wiki?name=SpatialIndex.
But neither the GeoPackage standard (nor its GDAL implementation) includes it.

Regarding the ROWID attribute, it is documented in https://www.sqlite.org/lang_createtable.html#rowid. In the GeoPackage standard, the Feature ID attribute refers to it.

The complete code I tested is the following:

import geopandas as gpd
import shapely.geometry
import numpy as np
import time

np.random.seed(42)

elapsed_with_rtree = []
elapsed_without_rtree = []

ns = [8,16, 24, 32, 40, 48, 56, 64]

for n in ns:

    p = gpd.GeoSeries(
        [shapely.geometry.box(j, i, j + 1, i + 1) for i in range(n) for j in range(n)]
    )

    q = gpd.GeoSeries(
        [
            shapely.geometry.Point(e)
            for e in np.random.uniform(low=0, high=n, size=[n * n * 10, 2])
        ]
    ).buffer(distance=0.1)
    
    !rm -f united.gpkg
    p.to_file("united.gpkg",layer="p")
    q.to_file("united.gpkg",layer="q")
    
    query = """
        SELECT ST_Difference(
          p.geom,
          (SELECT ST_UNION(geom)
           FROM q
           WHERE ST_Intersects(p.geom,geom)
           )
        ) AS geom FROM p;
    """
    query = ' '.join(query.split())
    
    rtree_query = """
        SELECT ST_Difference(
          p.geom,
          (SELECT ST_UNION(geom)
           FROM q
           WHERE ROWID IN(
             SELECT id
             FROM rtree_q_geom
             WHERE minx <= MbrMaxX(p.geom)
               AND maxx >= MbrMinX(p.geom)
               AND miny <= MbrMaxY(p.geom)
               AND maxy >= MbrMinY(p.geom)
             )
           AND ST_Intersects(p.geom, geom)
           )
        ) AS geom FROM p;
    """
    rtree_query = ' '.join(rtree_query.split())
    
    t0 = time.time()
    subprocess.run(['ogrinfo', '-q', '-sql', query, 'united.gpkg'], stdout=subprocess.DEVNULL)
    t1 = time.time()
    subprocess.run(['ogrinfo', '-q', '-sql', rtree_query, 'united.gpkg'], stdout=subprocess.DEVNULL)
    t2 = time.time()

    elapsed_without_rtree.append(t1-t0)
    elapsed_with_rtree.append(t2-t1)


import matplotlib.pyplot as plt

plt.figure(figsize=(8,4))

plt.scatter([e**2 for e in ns], elapsed_without_rtree, c='k')
plt.plot([e**2 for e in ns], elapsed_without_rtree, c='r', label="without rtree")

plt.scatter([e**2 for e in ns], elapsed_with_rtree, c='k')
plt.plot([e**2 for e in ns], elapsed_with_rtree, c='b', label="with rtree")

plt.xlabel("Number of geometries in p")
plt.ylabel("ogrinfo execution time")

plt.legend()