setting and getting bone transformations
Edit bone transforms are defined with head, tail and roll

Edit bones store the data in head, tail and roll. Applying the transform matrix as above to the edit bone sets the roll to 90 degrees. It does not translate the head or tail. Your choice of bone head and tail align to the z axis and the transform matrix is rotating around the z axis 90 degrees.

Note: To get the adjusted value of bone matrix need to toggle back into object mode

bone.transform(tr)
bpy.ops.object.mode_set() # back to object mode
print('=====================')


Running the code above without the transform M1 and with M2 and getting the bone matrix after toggling back into object mode. ,  calculate the transform from M1 to M2

&gt;&gt;&gt; M1 = Matrix([[1, 0, 0], [0, 0, -1], [0, 1, 0]])
&gt;&gt;&gt; M2 = Matrix([[0, 0, 1], [1, 0, 0], [0, 1, 0]])
&gt;&gt;&gt; T = M2 * M1.inverted()
&gt;&gt;&gt; T
Matrix(((0.0, -1.0, 0.0),
        (1.0, 0.0, 0.0),
        (0.0, 0.0, 1.0)))

&gt;&gt;&gt; T.to_euler()
Euler((0.0, -0.0, 1.5707963705062866), 'XYZ')


Remember that part of this matrix is the arbitrary choice of z axis for bone direction. In bone space the y axis is along the direction of the bone.  Notice how both bone matrices above will transform Y axis to Z. 

&gt;&gt;&gt; y = Vector((0, 1, 0))
&gt;&gt;&gt; M1 * y
Vector((0.0, 0.0, 1.0))

&gt;&gt;&gt; M2 * y
Vector((0.0, 0.0, 1.0))


To sum up, to calculate what transform has been applied to get from A to B, need A. A is not simply a "zero" in this case since you have defined head and tail before the transform.

If you make the bone tail above scalar * Vector((0, 1, 0)) ie Y axis direction, the bone matrix will be identity

Edit re comment

You don't get the "wrong" rotation you get an equivalent rotation. If you use the resulting bone matrix from above as your transform matrix, you will get same result for both. (the left hand rule.. or right can't remember) This goes back to my first comment under question. 

Eg you can rotate default cone to make it point from up to forward in a number of ways. -90 around x axis or 90 around y 90 around z, different transforms same result. In object mode we can see the transform. In edit mode, just like with the edit bone, only the result is saved, what transforms were used to create that result, even given the starting point this is unknown. 

Blender creates the bone matrix based on head and tail and roll. Trying to match back  the transform used from the result isn't 1:1

Let's investigate the two matrices from latest example

&gt;&gt;&gt; tr
Matrix(((0.0, 1.0, 0.0, 10.0),
        (1.0, 0.0, 0.0, 21.0),
        (0.0, 0.0, -1.0, 0.0),
        (-0.0, 0.0, -0.0, 1.0)))

&gt;&gt;&gt; tr2
Matrix(((0.0, 1.0, 0.0, 10.0),
        (-1.0, 0.0, 0.0, 21.0),
        (0.0, 0.0, 1.0, 0.0),
        (-0.0, 0.0, -0.0, 1.0)))

&gt;&gt;&gt; q = tr.to_quaternion()
&gt;&gt;&gt; q2 = tr2.to_quaternion()
&gt;&gt;&gt; q.rotation_difference(q2).to_euler()
Euler((-0.0, -3.1415927410125732, -0.0), 'XYZ')


Note the difference in rotation between the two is a full revolution in one axis,

Using the resultant bone matrix instead as transform matrix. 

===============
&lt;Matrix 4x4 ( 0.0000, 1.0000,  0.0000, 10.0000)
            (-1.0000, 0.0000,  0.0000, 21.0000)
            ( 0.0000, 0.0000,  1.0000,  0.0000)
            (-0.0000, 0.0000, -0.0000,  1.0000)&gt;
&lt;Matrix 4x4 ( 0.0000, 1.0000, 0.0000, 10.0000)
            (-1.0000, 0.0000, 0.0000, 21.0000)
            ( 0.0000, 0.0000, 1.0000,  0.0000)
            ( 0.0000, 0.0000, 0.0000,  1.0000)&gt;


How do you export bones relative to parent?

Exporting trouble. Bone tail and head positions change quaternion rotation