Instead of using the vector layer to calculate the zonal statistics of the raster (terra::extract() function), you can use raster data with zones (terra::zonal() function), which will be faster. Alternatively, you can use the very fast {exactextractr} package.

Code

{terra}

library("terra")
set.seed(1)
f = system.file("ex/lux.shp", package = "terra")
v = vect(f)
n = 3000
r = rast(v, nrows = n, ncols = n, names = "value", vals = rnorm(n ^ 2))
plot(r)
plot(v, add = TRUE)

system.time({
  zonal_vector = extract(r, v, fun = "mean")
})
##    user  system elapsed 
##    8.19    0.27    8.47
system.time({
  rr = rasterize(v, r, field = seq_len(nrow(v)))
  zonal_raster = zonal(r, rr, fun = "mean")
})
##    user  system elapsed 
##    0.18    0.02    0.18
all.equal(zonal_vector$value, zonal_raster$value)
## [1] TRUE
head(zonal_vector)

{exactextractr}

library("sf")
library("exactextractr")
v_sf = st_as_sf(v)
system.time({
  zonal_exactextractr = exact_extract(r, v_sf, fun = "mean", progress = FALSE)
})
##    user  system elapsed 
##    0.24    0.09    0.33
all.equal(zonal_exactextractr, zonal_vector$value, tolerance = 0.001)
## [1] TRUE