Find $\operatorname{trace}(BY^{-1})$, given $\mathrm{vec}(Y)=(A\otimes A-I_n)^{-1}\mathrm{vec}(B)$

Question

Given diagonal $A\in\mathbb{R}^{n\times n}$, symmetric $B\in\mathbb{R}^{n\times n}$, and $\mathrm{vec}(Y)=(A\otimes A-I_n)^{-1}\mathrm{vec}(B)$, find the following:

\begin{align} \operatorname{trace}(BY^{-1}). \end{align}

$A-I>0,$ $B$ is nonsingular.

---

My attempt:

For special case $A=aI_n$:

\begin{align} \mathrm{vec}(Y)=(A\otimes A-I_n)^{-1}\mathrm{vec}(B)=(aI_n\otimes aI_n-I_n)^{-1}\mathrm{vec}(B)=(a^2-1)^{-1}\mathrm{vec}(B). \end{align}

Thus

\begin{align} Y^{-1}=(a^2-1)B^{-1}, \end{align}

and \begin{align} \operatorname{trace}(BY^{-1})=\operatorname{trace}((a^2-1)I_n)=n(a^2-1). \end{align}

For special case $A=\begin{bmatrix} a_1 & & \\ & a_2 & \\ &&a_2 \end{bmatrix}:$

\begin{align} \operatorname{trace}(BY^{-1})=a_1^2a_2^2 + a_2^2 - 2. \end{align}

---

A=[3 0 0; 0 2 0; 0 0 2];
a=[3 2 2]';
T=zeros(1000,1);
for i=1:1000
    B=rand(3,2);
    [X,K,L] = idare(A,B,zeros(3,3));
    T(i)=a'*X.*inv(X)*a;
end
plot(T)

Answer 1

Define the vectors $$\eqalign{ a = {\rm diag}(A),\quad b = {\rm vec}(B),\quad y = {\rm vec}(Y) \\ }$$ and the matrices $$\eqalign{ A &= {\rm Diag}(a) \\ M &= A\otimes A - I_n\otimes I_n }$$ Use these to rearrange the given relationship into a matrix equation.
$$\eqalign{ y &= M^{-1}b \\ b &= My = \left(A\otimes A - I_n\otimes I_n\right)y &= {\rm vec}(AYA - Y) \\ B &= AYA - Y \\ BY^{-1} &= AYAY^{-1} - I_n \\ }$$ The function in question is the trace of that last expression. $$\eqalign{ \phi &= {\rm Tr}(BY^{-1}) \\&= {\rm Tr}(AYAY^{-1}) - {\rm Tr}(I_n) \\ &= YAY^{-1}:A - n \\ &= {\rm diag}(YAY^{-1}):{\rm diag}(A) - n \\ &= \left(Y^{-T}\odot Y\right)a:a - n \\ &= R_Ya:a - n \\ &= a^TR_Ya - n \\ }$$ The matrix $R_Y$ is known as the Relative Gain Array of $Y$ and has some interesting properties.
For example, denoting the all-ones vector by ${\tt1},\,$ it demonstrates Markov-like behavior
$$R_Y{\tt1}={\tt1},\qquad {\tt1}^TR_Y{\tt1}={\tt1}^T{\tt1}=n$$ Re-examining the special case of $\;\Big(A=\alpha I\iff a=\alpha{\tt1}\Big)\;$ yields $$\eqalign{ \phi &= (\alpha{\tt1})^TR_Y(\alpha{\tt1}) - n \\ &= \alpha^2n - n \\ \\ }$$

---

In the above, the symbol $(\odot)$ denotes the elementwise/Hadamard product
while a colon $(:)$ represents the trace/Frobenius product, i.e. $$\eqalign{ A:B = {\rm Tr}(A^TB) }$$ Finally, the diag() function creates a column vector from the main diagonal of the input matrix, while Diag() creates a diagonal matrix from a vector argument.

---

Update

I don't think you're interpreting your simulation correctly.

Here is a short snippet Julia/Matlab code which demonstrates that changing $B$ (while holding $A$ constant) definitely affects the value of $\phi$ and $R_Y$

a  = 100*rand(3,1);  A = Diag_(a)
   77.2603   0.0      0.0   
    0.0     83.9703   0.0   
    0.0      0.0     58.4881
     
B  = rand(3,3);  B += B';
y  = inv(kron(A,A)-I)*vec(B);  Y = reshape(y,size(B));
Ry = inv(Y').*Y
   0.564956   2.99807   -2.56303
   2.99807   -0.343684  -1.65439
  -2.56303   -1.65439    5.21742

a'*Ry*a
   18283.48744321049

B  = rand(3,3);  B += B';
y  = inv(kron(A,A)-I)*vec(B);  Y = reshape(y,size(B));
Ry = inv(Y').*Y
3×3 Array{Float64,2}:
   21.8927   17.0631   -37.9558
   17.0631    1.88734  -17.9505
  -37.9558  -17.9505    56.9063

a'*Ry*a
   40704.04750730478