This paper investigates the joint design of multiple non-regenerative multiple-input multiple-output (MIMO) relaying matrices, with the purpose of minimizing the mean square error (MSE) between the transmitted signals from the source and the received signals at the destination. Two types of constraints on the transmit power of the relays are considered separately: 1) a weighted sum power constraint, and 2) per-relay power constraints. As opposed to using general-purpose interior-point methods, we exploit the inherent structure of the problems to develop more efficient algorithms. Under the weighted sum power constraint, the optimal solution is expressed as a function of a Lagrangian parameter. By introducing a complex scaling factor at the destination, we derive a closed-form expression for this parameter, thereby avoiding the need to solve an implicit nonlinear equation numerically. Under the per-relay power constraints, the optimal solution is the same as that under the weighted sum power constraint if particular weights are chosen. We then propose an iterative power balancing algorithm to compute these weights. In addition, under both types of constraints, we investigate the joint design of a MIMO equalizer at the destination and the relaying matrices, using block coordinate descent or steepest descent. The bit-error rate (BER) simulation results demonstrate that all the proposed designs, under either type of constraints, with or without the equalizer, perform much better than previous methods.