Qiskit provides the qiskit.quantum_info.Operator class to get the unitary matrix operator from the corresponding quantum circuit, as in the following example:

from qiskit import QuantumCircuit
from qiskit.quantum_info import Operator
from qiskit.visualization import array_to_latex

qc = QuantumCircuit(2)
qc.h(0)

op = Operator(qc)
array_to_latex(op)

\begin{bmatrix} \frac{1}{\sqrt2} & \frac{1}{\sqrt2} & 0 & 0\\ \frac{1}{\sqrt2} & -\frac{1}{\sqrt2} & 0 & 0\\ 0 & 0 & \frac{1}{\sqrt2} & \frac{1}{\sqrt2} \\ 0 & 0 & \frac{1}{\sqrt2} & -\frac{1}{\sqrt2} \end{bmatrix}

However, using Operator(QuantumCircuit) throws an error in the case of a parameterized quantum circuit (PQC):

from qiskit.circuit import Parameter

qc = QuantumCircuit(2)
theta = Parameter(name='$\\theta$')
qc.ry(theta, 0)

op = Operator(qc)  # ERROR!

This brings me to the question: is there a way in Qiskit to get the matrix operator symbolic representation from a given arbitrary PQC? For instance, in this case I would like to get a sympy.matrices.dense.Matrix object (with just one parameter $\theta$) like this:

\begin{bmatrix} \cos\left(\frac{\theta}{2}\right) & -\sin\left(\frac{\theta}{2}\right) & 0 & 0\\ \sin\left(\frac{\theta}{2}\right) & \cos\left(\frac{\theta}{2}\right) & 0 & 0\\ 0 & 0 & \cos\left(\frac{\theta}{2}\right) & -\sin\left(\frac{\theta}{2}\right)\\ 0 & 0 & \sin\left(\frac{\theta}{2}\right) & \cos\left(\frac{\theta}{2}\right) \end{bmatrix}