Thank you, but this was not really what I asked. Anyway, I figured out the answer myself: in the axis-angle representation, the rotation axis is perpendicular to both the world bone axis and the y-axis (= local bone axis), and the angle is given by cos(phi) = dot product. Here is a code excerpt, using transformations.py.

```
def matrixLocalFromBone(self):
u = self.tail.sub(self.head)
length = sqrt(u.dot(u))
if length < 1e-3:
print("Zero-length bone %s. Removed" % self.name)
self.matrix_local = tm.identity(4)
self.matrix_local.matrix[:3,3] = self.head.vector
return
u = u.div(length)
ey = Vector((0,1,0))
yu = ey.dot(u)
if abs(yu) > 0.999:
axis = ey
if yu > 0:
angle = 0
else:
angle = pi
else:
axis = ey.cross(u)
length = sqrt(axis.dot(axis))
axis = axis.div(length)
angle = acos(yu)
mat = tm.rotation_matrix(angle,axis)
matrix = Matrix(mat)
if self.parent:
pmat = self.parent.matrix_local.inverted()
self.matrix_local = matrix.mult(pmat)
else:
self.matrix_local = matrix
self.matrix_local.matrix[:3,3] = self.head.vector
```

Actually, this is the answer for roll = 0. For nonzero roll, one must multiply with a extra rotation around the y axis somewhere.