Q:

How is the precision of a quantum simulation algorithm actually proved? 

The problem of quantum simulation can be formulated as follows:

Given a Hamiltonian $H$ ( $2^n \times 2^n$ hermitian matrix acting on $n$ qubits), a time $t$ and maximum simulation error $\epsilon$, the goal is to find an algorithm that approximates $U$ such that $ || U - e^{-iHt}|| \leq \epsilon $, where $ e^{-iHt}$ is the ideal evolution and $|| \cdot ||$ is the spectral norm. 

I was wondering, how can you compare the result of your algorithm with an ideal evolution?  It would make sense, if this ideal evolution could be calculated analytically, but I think it's rarely the case.  This $e^{-iHt}$ could be also given by an oracle, but in a real world scenario we just cannot use such an approach.  
So, how can you actually compare the result of your approach with something, that cannot be efficiently obtained (in fact, if it could be obtained efficiently, then the whole idea of designing other algorithm to approximate this result doesn't make sense anymore). 

A:

how can you compare the result of your algorithm with an ideal evolution? 

You cannot and you do not need to. 
As you said, computing $e^{-iHt}$ is intractable for most of the interesting cases.  If it was not, chemistry simulations would be easy, solving the Schrödinger equation too. 
The thing you can do though is to prove that your algorithm will, for a (potentially restricted) set of Hamiltonians, output a unitary $U$ that satisfy the condition $||U - e^{-iHt}|| \leqslant \epsilon$. You do not compare the output of your algorithm with the ideal evolution, you prove that your algorithm is guaranteed to return a result that satisfy the condition. 
If you are more familiar with computer science, these are pre/post-conditions.  
Applied to the definition of Hamiltonian simulation you gave:

Pre-conditions are:

The matrix given should be Hermitian. 
The matrix given should be of size $\left( 2^n, 2^n \right)$ with $n \in \mathbb{N}^*$.
$t \in \mathbb{R}$ (note that I did not excluded negative times, it might be an error from my side). 
$\epsilon \in \mathbb{R}^*_+$.

Post-conditions are:

$U$ is a quantum circuit composed of quantum gates. 
The unitary matrix "implemented" by $U$ satisfy $||U - e^{-iHt}|| \leqslant \epsilon$. 

Any algorithm that solve the Hamiltonian Simulation problem should, when the pre-conditions are checked (i.e. the user gave correct inputs), return a result that verify the post-conditions above.