Computing with matrix differential operators

Here is a modification of your code that works. The problem was that you weren't combining the Compositions for each matrix entry correctly in CircleTimes: using Plus with Composition doesn't lead to a function which can be applied subsequently to another expression. You have to make each sum of compositions into a function itself.

First I'll copy the definitions that I didn't modify:

Clear[f, g, x];

I1 := {{0 &, (x # + D[#, x] &)}, {0 &, -D[#, x] &} }

I2 := {{0 &, x^2 # &}, {0 &, D[#, x] &}}

u = {f[x], g[x]};

CenterDot[a_, v_] := Inner[Through[#1[#2] &, Plus], a, v]

Now my modifications:

Clear[CircleTimes];
CircleTimes[a_, b_] := 
 Inner[Composition[#1, #2] &, a, b, 
  Function[f, Total@Through[{##}[f]]] &]


CirclePlus[a_, b_] := 
 MapThread[Function[f, Total@Through[{##}[f]]] &, {a, b}, 2]

The modification in CircleTimes is that I specified the replacement of the Plus operation in Inner as Function[f, Total@Through[{##}[f]]] & which is now itself a function that adds the individual operators after applying them to whatever argument f you provide.

Next, the same issue with Plus also arises in I3. So I defined an addition of operators CirclePlus in which the same trick is used as above, by combining the matrices of operators element-wise using MapThread.

Now the tests:

CenterDot[I1, u]

(* ==> {x g[x] + Derivative[1][g][x], -Derivative[1][g][x]} *)

CenterDot[I2, u]

(* ==> {x^2 g[x], Derivative[1][g][x]} *)

CenterDot[CircleTimes[I1, I2], u]

(*
==> {x Derivative[1][g][x] + (g^\[Prime]\[Prime])[x], -(
   g^\[Prime]\[Prime])[x]}
*)

I3 := CirclePlus[CircleTimes[I1, I2], CircleTimes[I2, I1]]

CenterDot[I3, u]

(*
==> {x Derivative[1][g][x] - 
  x^2 Derivative[1][g][x] + (g^\[Prime]\[Prime])[x], -2 (
   g^\[Prime]\[Prime])[x]}
*)

Id := {{1 &, 0 &}, {0 &, 1 &}}

CenterDot[CirclePlus[Id, Id], u]

(* ==> {2, 2} *)

As you can see, the tests produce outcomes that look correct. Edit: the last example is from the comment, and I changed the code to accommodate such constant operators. Note that this is not the same as the identity which multiplies the test function by 1. My understanding of the comment was that the desired operator was meant to replace the test function by 1.

This is how I would define your functions, similar but I think simpler than @Jens answer:

ClearAll[CirclePlus, CircleTimes, CenterDot]

SetAttributes[CirclePlus, {Flat, Listable}]

CirclePlus[a__][x_] := Activate @* Through @* Inactive[Plus][a] @ x
CircleTimes[a_, b_] := Inner[Composition, a, b, CirclePlus]
CenterDot[a_, b_] := Inner[Compose, a, b]

For the OP example:

I1 := {{0&, (x # + D[#,x]&)}, {0&, -D[#,x]&}}
I2 := {{0&, x^2 #&}, {0&, D[#,x]&}}
u = {f[x], g[x]};

The same results as @Jens:

CenterDot[I1, u] //TeXForm

$\left\{g'(x)+x g(x),-g'(x)\right\}$

CenterDot[I2, u] //TeXForm

$\left\{x^2 g(x),g'(x)\right\}$

CenterDot[CircleTimes[I1, I2], u] //TeXForm

$\left\{g''(x)+x g'(x),-g''(x)\right\}$

I3 := CirclePlus[CircleTimes[I1, I2], CircleTimes[I2, I1]]
CenterDot[I3, u] //TeXForm

$\left\{g''(x)+x^2 \left(-g'(x)\right)+x g'(x),-2 g''(x)\right\}$

Id := {{1&, 0&}, {0&, 1&}}
CenterDot[CirclePlus[Id, Id], u] //TeXForm

$\{2,2\}$