how about we use same method for invert(A) as for A^^-1 ?

Another example:

  (%i1) (matrix_element_mult : ".", matrix_element_transpose : 'transpose)$

  (%i2) m : matrix([a,b],[c,d])$

  (%i3) mi : invert_by_lu(m,'noncommutingring);

  (%o3) matrix([a^^(-1) . (b . (d-c . a^^(-1) . b)^^(-1) . c . a^^(-1)+1),-a^^(-1) . b . (d-c . a^^(-1) . b)^^(-1)],[-(d-c . a^^(-1) . b)^^(-1) . c . a^^(-1),(d-c . a^^(-1) . b)^^(-1)])

I'm not sure why ?subst does what I want (doesn't  substitute  for -1, but it does substitute for 1), but I think it's OK

  (%i4) mi : ?subst(I,1,%);
  (%o4) matrix([a^^(-1) . (b . (d-c . a^^(-1) . b)^^(-1) . c . a^^(-1)+I),-a^^(-1) . b . (d-c . a^^(-1) . b)^^(-1)],[-(d-c . a^^(-1) . b)^^(-1) . c . a^^(-1),(d-c . a^^(-1) . b)^^(-1)])

Let's subst 2x2 matrices for a thru d and I

  (%i5) subs : [I = matrix([1,0],[0,1]), a = matrix([3,4],[5,7]), b = matrix([6,7],[9,21]), c = matrix([1,1],[0,1]), d = matrix([9,1],[5,6])]$

Either use psubst or declare a,b,c,d,I to be nonscalar--use subst and you'll get a big  mess (but maybe correct--I don't know. Try to predict the shape of the
result before trying it):

  (%i7) mii : psubst(subs,mi);
  (%o7) matrix([matrix([-107/98,69/98],[74/49,-73/98]),matrix([37/49,-129/98],[-79/98,48/49])],[matrix([-1/49,-1/98],[-23/98,13/98]),matrix([11/98,2/49],[2/49,-3/98])])

Checks:

  (%i8) subst(subs,m) . mii;
  (%o8) matrix([matrix([1,0],[0,1]),matrix([0,0],[0,0])],[matrix([0,0],[0,0]),matrix([1,0],[0,1])])

  (%i9) mii . subst(subs,m);
  (%o9) matrix([matrix([1,0],[0,1]),matrix([0,0],[0,0])],[matrix([0,0],[0,0]),matrix([1,0],[0,1])])

--Barton