Without further assumptions, I think $F$ is not necessarily concave.
Let $\mathbf{X}_1\sim p_1$, $\mathbf{X}_2\sim p_2$ and $B\sim\textrm{Bernoulli}(\lambda)$ be independent, and let
\begin{align*}
\mathbf{X} &:=
\begin{cases}
\mathbf{X}_1 & \text{if $B=1$,} \\
\mathbf{X}_2 & \text{if $B=0$.}
\end{cases}
\end{align*}
Then, $\mathbf{X}\sim p=\lambda p_1 + (1-\lambda) p_2$.
In general, for two random variables $Z$ and $C$, where $Z$ is continuous and $C$ is discrete, we have
\begin{align*}
h(Z) + H(C\,|\,Z) &= H(C) + h(Z\,|\,C) \;,
\end{align*}
where $H(\cdot)$ denotes the ordinary (discrete) entropy and $h(\cdot)$ is the differential entropy.
It follows that
\begin{align*}
&
\overbrace{h(\mathrm{A}\mathbf{X}) - h(\mathrm{A}\mathrm{U}\mathbf{X})}^{F(p)}
+
\overbrace{H(B\,|\,\mathrm{A}\mathbf{X})
- H(B\,|\,\mathrm{A}\mathrm{U}\mathbf{X})}^{\displaystyle(\sharp)} \\
&=
h(\mathrm{A}\mathbf{X}\,|\,B)
- h(\mathrm{A}\mathrm{U}\mathbf{X}\,|\,B)
+ H(B) - H(B) \\
&=
\lambda\big(\underbrace{h(\mathrm{A}\mathbf{X}_1) - h(\mathrm{A}\mathrm{U}\mathbf{X}_1)}_{F(p_1)}\big)
+
(1-\lambda)\big(\underbrace{h(\mathrm{A}\mathbf{X}_2) - h(\mathrm{A}\mathrm{U}\mathbf{X}_2)}_{F(p_2)}\big)
\end{align*}
provided that $p_1$ and $p_2$ are absolutely continuous w.r.t. the three-dimensional Lebesgue and $\mathrm{A}$ is almost surely non-singular. (Otherwise, the differential entropies become $-\infty$ and $F$ would not be well-defined.)
Therefore, in order for $F$ to be concave, we must have
\begin{align*}
(\sharp) = H(B\,|\,\mathrm{A}\mathbf{X})
- H(B\,|\,\mathrm{A}\mathrm{U}\mathbf{X})
&\leq 0 \tag{?}
\end{align*}
whenever $p_1$ and $p_2$ are absolutely continuous and $\mathrm{A}$ is almost surely non-singular.
[Update: The original example was not valid because it disregarded the requirement that $p_1$ and $p_2$ have to be absolutely continuous and $\mathrm{A}$ non-singular. The following sketch is meant to circumvent that issue.]
Fix $0<\lambda<1$. Let
\begin{align*}
\hat{\mathrm{A}} &:=
\begin{bmatrix}
1 & 1/2 & 1/2 \\
0 & -1/2 & 1/2 \\
0 & -1/2 & 1/2
\end{bmatrix}
&
\hat{\mathbf{X}}_1 &:=
\begin{bmatrix}
1 \\ 0 \\ 0
\end{bmatrix}
&
\hat{\mathbf{X}}_2 &:=
\begin{bmatrix}
0 \\ 1 \\ 1
\end{bmatrix}
\end{align*}
Let $\mathrm{A}$ be a non-singular (deterministic or random) matrix which is very close to $\hat{\mathrm{A}}$, and let $\mathbf{X}_1=\hat{\mathbf{X}}+\sigma\mathbf{Z}_1$ and $\mathbf{X}_2=\hat{\mathbf{X}}+\sigma\mathbf{Z}_2$, where $\mathbf{Z}_1$ and $\mathbf{Z}_2$ are two independent standard normal vectors and $\sigma$ is very small. Assume that $\mathbf{Z}_1$, $\mathbf{Z}_2$, $\mathrm{U}$ and $\mathrm{A}$ are all independent.
Note that both $\mathrm{A}\mathbf{X}_1$ and $\mathrm{A}\mathbf{X}_2$ are highly concentrated around a vector very close to $\hat{\mathbf{X}}_1$. By chooseing $\mathrm{A}$ close enough to $\hat{\mathrm{A}}$, we can make sure that $\mathrm{A}\mathbf{X}_1$ and $\mathrm{A}\mathbf{X}_2$ are hardly distinguishable. Hence, $\mathrm{A}\mathbf{X}$ would hardly have any information about $B$, and as a result
\begin{align*}
H(B\,|\,\mathrm{A}\mathbf{X}) &\approx H(B) = H(\lambda) \;.
\end{align*}
On the other hand, $\mathrm{A}\mathrm{U}\mathbf{X}_1$ and $\mathrm{A}\mathrm{U}\mathbf{X}_2$ will be distinguishable, with $\mathrm{A}\mathrm{U}\mathbf{X}_1$ still being close to the linear span of $\hat{\mathbf{X}}_1$ and $\mathrm{A}\mathrm{U}\mathbf{X}_2$ typically far from it. In particular, $\mathrm{A}\mathrm{U}\mathbf{X}$ has significant information about $B$ and hence
\begin{align*}
H(B\,|\,\mathrm{A}\mathrm{U}\mathbf{X}) &\ll H(B) = H(\lambda) \;.
\end{align*}
Therefore, in this example, $(\sharp)>0$ contrary to the claim.
