The resting yaw rotation is determined by the magnetometer - it's the only one of the three sensors that has an absolute reference direction in the horizontal plane.
Unfortunately, it's also the most fiddly of the sensors. Metal components inside the device can cause its detected magnetic field vector to be shifted (hard iron distortion) or stretched/rotated/skewed (soft iron distortion).
The most robust way to solve this is to calibrate the device prior to use.
If you ask the user to tumble the sensor around in all directions, the magnetic field vector that you read from the magnetometer will trace out an ellipsoid. Once you have enough samples around the whole surface of the ellipsoid, you can estimate its shape. You can then compute a minimal transformation that maps the ellipsoid to a unit sphere centered at the origin.
This can be as simple as keeping track of the min & max you've observed on each axis, shifting the result by
- (min+max)/2 and scaling it by
2/(max - min), though if you have substantial soft iron distortions that skew/rotate your ellipsoid, or misalignments of the sensor chip itself, you may need a more sophisticated correction.
Now you can apply this calibrating transformation to your magnetic field samples before passing them as input to your sensor-fusion routine to estimate the device orientation. This should minimize the differences in the two devices' estimates.
A simpler but less robust calibration you can try is to ask the user to set the two devices on a flat surface in a standard orientation (it helps if the devices have a flat bottom to rest on, and a clear line somewhere to point forward).
Wait a moment for the orientation estimate to settle, then compare the quaternions you get from each device.
// Total difference in orientation.
// (Optionally, you can filter this to consider only the difference in yaw)
Quaternion AtoB = deviceB.rotation * Quaternion.Inverse(deviceA.rotation);
// We'll average-out the two and meet in the middle:
// correcting A halfway to B's orientation,
// and B halfway to A's.
Quaternion aCorrection = Quaternion.Lerp(Quaternion.identity, AtoB, 0.5f);
Quaternion bCorrection = Quaternion.Inverse(aCorrection);
Now you can apply this correction as
Quaternion correctedA = aCorrection * deviceA.orientation;
Quaternion correctedB = bCorrection * deviceB.orientation;
The reason I say this is less robust is because it assumes the same yaw offset applies all the time - but this can be sensitive to the orientation of the device. You could repeat the calibration process in multiple orientations (standing on end, on its side, etc...), saving multiple different correction quaternions, and then blend between them depending on the current device orientation. We're still treating the symptom here (differing yaw estimates) rather than the root cause (differently-biased magnetometer readings).