If I from qiskit.quantum_info import random_unitary and then random_unitary(2**number_of_qubits) What returns is a unitary matrix with dimension (2^number_of_qubits,2^number_of_qubits), and if I stimulating the unitary operation on a classical computer it is matrix computation.

Here comes my question: what quantum computers do to implement this random unitary?

From the universality of quantum gate set, we know that an arbitrary single-qubit operation can be decomposed into the set of {H, S, T} or $\hat U=\hat\Phi(\delta)\hat R_z(\alpha)\hat R_y(\theta)\hat R_z(\beta)$ where $\hat \Phi(\delta)=\begin{pmatrix}e^{i\delta}&0\\0&e^{i\delta}\end{pmatrix}$, $\hat R_z(\alpha)=\begin{pmatrix}e^{i\alpha/2}&0\\0&e^{-i\alpha/2}\end{pmatrix}$ and $\hat R_y(\theta)=\begin{pmatrix}cos\frac{\theta}{2}&sin\frac{\theta}{2}\\-sin\frac{\theta}{2}&cos\frac{\theta}{2}\end{pmatrix}$.

But when considering physical realizations, we need to consider the problem of finite precision and so on. Then, which method will the quantum computers adopt? And more generally, what if we introduce CNOT here to form the universal quantum gate set for an arbitrary number of qubits?

Even a great reference paper or detailed qiskit documentation can be helpful(in fact to make my work rigorous this is essential). I appreciate you in advance.

If I from qiskit.quantum_info import random_unitary and then random_unitary(2**number_of_qubits) What returns is a unitary matrix with dimension $(2^{number\_of\_qubits},2^{number\_of\_qubits})$, and if I stimulating the unitary operation on a classical computer it is matrix computation.

Here comes my question: what quantum computers do to implement this random unitary?

From the universality of quantum gate set, we know that an arbitrary single-qubit operation can be decomposed into the set of $\{H, S, T\}$ or $\hat U=\hat\Phi(\delta)\hat R_z(\alpha)\hat R_y(\theta)\hat R_z(\beta)$ where $\hat \Phi(\delta)=\begin{pmatrix}e^{i\delta}&0\\0&e^{i\delta}\end{pmatrix}$, $\hat R_z(\alpha)=\begin{pmatrix}e^{i\alpha/2}&0\\0&e^{-i\alpha/2}\end{pmatrix}$ and $\hat R_y(\theta)=\begin{pmatrix}cos\frac{\theta}{2}&sin\frac{\theta}{2}\\-sin\frac{\theta}{2}&cos\frac{\theta}{2}\end{pmatrix}$.

But when considering physical realizations, we need to consider the problem of finite precision and so on. Then, which method will the quantum computers adopt? And more generally, what if we introduce CNOT here to form the universal quantum gate set for an arbitrary number of qubits?

Even a great reference paper or detailed qiskit documentation can be helpful  (in fact to make my work rigorous this is essential). I appreciate you in advance.