Can you work with homogeneous coordinates in the Python API? If I do this (excuse the transpose, I'm used to looking at matrices the other way round):

import bpy
from  mathutils import Matrix
ob =  bpy.context.object

m = Matrix( ([1, 0, 0, 0],
             [0, 1, 0, 0],
             [0, 0, 1, 0],
             [0, 0, 1, 0] )  )

m.transpose()

for vert in ob.data.vertices:
    vert.co = m @ vert.co


I get a translation of the mesh by 1 in object-Z. Good. So somewhere, vert.co must be understood to be (x,y,z,w=1), otherwise the translation wouldn't be picked up.

However, if I do this:

import bpy
from  mathutils import Matrix
ob =  bpy.context.object

m = Matrix( ([1, 0, 0, 0],
             [0, 1, 0, 0],
             [0, 0, 1, 1],
             [0, 0, 0, 0] )  )

m.transpose()

for vert in ob.data.vertices:
    vert.co = m @ vert.co


.. I'm hoping for a basic perspective transform, in which the w of every coordinate is set to its z, and then, in normalization, the whole coordinate is divided through by its w before using its x,y,z. But, OK, nothing happens.

Is there any way of getting at the implicit 'w' of the homogeneous vector, or do I have to divide it by hand? 

Is there a way using matrix multiplication like this to make projections in the API? 

It could be my Blender, it could be my math. I don't know....
Make the vert.co 4d


  From Perspective projection of transformation matrix
  wiki
  
  The simplest perspective projection uses the origin as the center of
  projection, and the plane at ${\displaystyle z=1}$  as the image
  plane. The functional form of this transformation is then 
  ${\displaystyle x'=x/z}$ ; ${\displaystyle y'=y/z}$. We can express
  this in homogeneous coordinates as: 
  
  $\begin{pmatrix}x_c\\ y_c \\ z_c \\ w_c\end{pmatrix} = \begin{pmatrix}1 &amp; 0 &amp; 0 &amp; 0\\ 0 &amp; 1 &amp; 0 &amp; 0 \\ 0 &amp; 0 &amp; 1 &amp; 0 \\ 0 &amp; 0 &amp; 1 &amp; 0\end{pmatrix} \begin{pmatrix}x\\ y \\ z \\ w=1\end{pmatrix}$
  
  $\begin{pmatrix}x'\\ y' \\ z' \\ 1\end{pmatrix} = 1 / w_c \begin{pmatrix}x_c\\ y_c \\ z_c \\ w_c\end{pmatrix}$


To emulate above would map the vert co to 4d using w of the result


  Is there any way of getting at the implicit 'w' of the homogeneous
  vector, or do I have to divide it by hand?


yeah suppose I've done that.

import bpy
from  mathutils import Matrix
from math import radians
ob =  bpy.context.object

m = Matrix( ([1, 0, 0, 0],
             [0, 1, 0, 0],
             [0, 0, 1, 1],
             [0, 0, 0, 0] )  )

m.transpose()

R = Matrix.Rotation(radians(30), 4, 'X')
#m = R @ m
for vert in ob.data.vertices:
    v = m @ vert.co.to_4d() # sets vert.co.w to 1

    '''
    # vert.co.w is set to vert.co.z
    v =  vert.co.to_4d()
    v.w = v.z
    v = m @ v
    '''
    vert.co = 1 / v.w * v.to_3d()


If run on default cube, flattens it to the plane z=1.  However others will fail with zero divide error.

The question is how to handle $\ w_c = 0$


  One property of homogeneous coordinates is that they allow you to have
  points at infinity (infinite length vectors), which is not possible
  with 3D coordinates. Points at infinity occur when W=0. If you try and
  convert a W=0 homogeneous coordinate into a normal W=1 coordinate, it
  results in a bunch of divide-by-zero operations:
  
  This means that homogeneous coordinates with W=0
  can not be converted into 3D coordinates.
  What use does this have? Well, directional lights can be though of as
  point lights that are infinitely far away. When a point light is
  infinitely far away the rays of light become parallel, and all of the
  light travels in a single direction. This is basically the definition
  of a directional light.


Playing around with meshes using above Not sure how useful the answer is. In first example in question above, worth noting m.translation = (0, 0, 1) whereas in second it is zero.  Converting 3d vector to 4d sets w to 1.  Which I assume is done for us since it doesn't spit the dummy multiplying a 4x4 with a 3x1.