C# PositiveDefiniteMatrix类代码示例

 /// <summary>
 /// VMP message to 'B'.
 /// </summary>
 /// <param name="X">Incoming message from 'X'. Must be a proper distribution.  If any element is uniform, the result will be uniform.</param>
 /// <param name="A">Incoming message from 'A'.</param>
 /// <param name="result">Modified to contain the outgoing message.</param>
 /// <returns><paramref name="result"/></returns>
 /// <remarks><para>
 /// The outgoing message is the exponential of the integral of the log-factor times incoming messages, over all arguments except 'A'.
 /// The formula is <c>int log(f(A,x)) q(x) dx</c> where <c>x = (X,B)</c>.
 /// </para></remarks>
 /// <exception cref="ImproperMessageException"><paramref name="X"/> is not a proper distribution</exception>
 /// <exception cref="ImproperMessageException"><paramref name="B"/> is not a proper distribution</exception>
 public static VectorGaussianArray BAverageLogarithm([SkipIfAllUniform] GaussianArray2D X, [SkipIfAllUniform] VectorGaussianArray A, VectorGaussianArray result)
 {
     int I = X.GetLength(0), J = X.GetLength(1), K = A[0].Dimension;
     var Eaap = new PositiveDefiniteMatrix[I];
     var ma = new Vector[I];
     for (int i = 0; i < I; i++)
     {
         Eaap[i] = new PositiveDefiniteMatrix(K, K);
         ma[i] = Vector.Zero(K);
         A[i].GetMeanAndVariance(ma[i], Eaap[i]);
         Eaap[i].SetToSumWithOuter(Eaap[i], 1, ma[i], ma[i]);
     }
     for (int j = 0; j < J; j++)
     {
         result[j].Precision.SetAllElementsTo(0);
         result[j].MeanTimesPrecision.SetAllElementsTo(0);
         for (int i = 0; i < I; i++)
         {
             // nb: would be more memory efficient to have a SetToAPlusCB routine
             result[j].Precision.SetToSum(result[j].Precision, Eaap[i] * X[i, j].Precision);
             result[j].MeanTimesPrecision.SetToSum(result[j].MeanTimesPrecision, ma[i] * X[i, j].MeanTimesPrecision);
         }
     }
     return result;
 }

		/// <summary>
		/// VMP message to 'innerProduct'
		/// </summary>
		/// <param name="A">Constant value for 'a'.</param>
		/// <param name="BMean">Buffer 'BMean'.</param>
		/// <param name="BVariance">Buffer 'BVariance'.</param>
		/// <returns>The outgoing VMP message to the 'innerProduct' argument</returns>
		/// <remarks><para>
		/// The outgoing message is the factor viewed as a function of 'innerProduct' conditioned on the given values.
		/// </para></remarks>
		public static Gaussian InnerProductAverageLogarithm(Vector A, Vector BMean, PositiveDefiniteMatrix BVariance)
		{
			Gaussian result = new Gaussian();
			// Uses John Winn's rule for deterministic factors.
			// Strict variational inference would set the variance to 0.
			// p(x) = N(a' E[b], a' var(b) a)
			result.SetMeanAndVariance(A.Inner(BMean), BVariance.QuadraticForm(A));
			return result;
		}

		/// <summary>
		/// Update the buffer 'CovarianceOfB'
		/// </summary>
		/// <param name="B">Incoming message from 'B'. Must be a proper distribution.  If any element is uniform, the result will be uniform.</param>
		/// <returns>New value of buffer 'CovarianceOfB'</returns>
		/// <exception cref="ImproperMessageException"><paramref name="B"/> is not a proper distribution</exception>
		public static PositiveDefiniteMatrix Eaat(GaussianArray A, PositiveDefiniteMatrix result)
		{
			int inner = A.Count;
			var VarianceOfA = Vector.Zero(inner);
			var MeanOfA = Vector.Zero(inner);
			for (int k = 0; k < inner; k++) {
				MeanOfA[k] = A[k].GetMean();
				VarianceOfA[k] = A[k].GetVariance();
			}
			result.SetToDiagonal(VarianceOfA);
			result.SetToSumWithOuter(result, 1, MeanOfA, MeanOfA);
			return result;
		}

        // Actually it doesn't require matrix to be positive semidefinite
        public static PositiveDefiniteMatrix Invert(PositiveDefiniteMatrix matrix)
        {
            Debug.Assert(matrix.Rows == 2 && matrix.Cols == 2);

            double determinant = matrix[0, 0] * matrix[1, 1] - matrix[0, 1] * matrix[1, 0];
            Debug.Assert(Math.Abs(determinant) > 1e-10);

            double invDeterminant = 1.0 / determinant;
            PositiveDefiniteMatrix result = new PositiveDefiniteMatrix(2, 2);
            result[0, 0] = invDeterminant * matrix[1, 1];
            result[0, 1] = -invDeterminant * matrix[0, 1];
            result[1, 0] = -invDeterminant * matrix[1, 0];
            result[1, 1] = invDeterminant * matrix[0, 0];

            return result;
        }

		public static Gamma BAverageConditional(PositiveDefiniteMatrix Product, PositiveDefiniteMatrix A)
		{
			if (Product.Count == 0) return Gamma.Uniform();
			bool allZeroA = true;
			double ratio = 0;
			for (int i = 0; i < Product.Count; i++) {
				if (A[i] != 0) {
					ratio = Product[i]/A[i];
					allZeroA = false;
				}
			}
			if (allZeroA) return Gamma.Uniform();
			for (int i = 0; i < Product.Count; i++) {
				if (Math.Abs(Product[i] - A[i]*ratio) > 1e-15) throw new ConstraintViolatedException("Product is not a multiple of B");
			}
			return Gamma.PointMass(ratio);
		}

        public static VectorGaussianWishart Combine(VectorGaussian position, Wishart orientation, VectorGaussianWishart result)
        {
            if (orientation.IsUniform())
            {
                result.SetToUniform();
            }
            else if (position.IsUniform())
            {
                result.SetTo(orientation.Shape, orientation.Rate, Vector.Zero(2), 0);
            }
            else
            {
                PositiveDefiniteMatrix rateTimesPrecision = new PositiveDefiniteMatrix(2, 2);
                rateTimesPrecision.SetToProduct(orientation.Rate, position.Precision);
                double trace = MathHelpers.Invert(rateTimesPrecision).Trace();
                Vector positionMean = position.MeanTimesPrecision * MathHelpers.Invert(position.Precision);
                result.SetTo(orientation.Shape, orientation.Rate, positionMean, orientation.Dimension / (orientation.Shape * trace));
            }

            return result;
        }

		/// <summary>
		/// Generates a data set from a particular true model.
		/// </summary>
		public Vector[] GenerateData(int nData)
		{
			Vector trueM1 = Vector.FromArray(2.0, 3.0);
			Vector trueM2 = Vector.FromArray(7.0, 5.0);
			PositiveDefiniteMatrix trueP1 = new PositiveDefiniteMatrix(
				new double[,] { { 3.0, 0.2 }, { 0.2, 2.0 } });
			PositiveDefiniteMatrix trueP2 = new PositiveDefiniteMatrix(
				new double[,] { { 2.0, 0.4 }, { 0.4, 4.0 } });
			VectorGaussian trueVG1 = VectorGaussian.FromMeanAndPrecision(trueM1, trueP1);
			VectorGaussian trueVG2 = VectorGaussian.FromMeanAndPrecision(trueM2, trueP2);
			double truePi = 0.6;
			Bernoulli trueB = new Bernoulli(truePi);
			// Restart the infer.NET random number generator
			Rand.Restart(12347);
			Vector[] data = new Vector[nData];
			for (int j = 0; j < nData; j++) {
				bool bSamp = trueB.Sample();
				data[j] = bSamp ? trueVG1.Sample() : trueVG2.Sample();
			}
			return data;
		}

		/// <summary>
		/// Update the buffer 'BVariance'
		/// </summary>
		/// <param name="B">Incoming message from 'b'. Must be a proper distribution.  If any element is uniform, the result will be uniform.</param>
		/// <param name="result">Modified to contain the outgoing message</param>
		/// <returns><paramref name="result"/></returns>
		/// <remarks><para>
		/// 
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="B"/> is not a proper distribution</exception>
		public static PositiveDefiniteMatrix BVariance([Proper] VectorGaussian B, PositiveDefiniteMatrix result)
		{
			return B.GetVariance(result);
		}

		/// <summary>
		/// EP message to 'product'
		/// </summary>
		/// <param name="A">Constant value for 'a'.</param>
		/// <param name="BMean">Buffer 'BMean'.</param>
		/// <param name="BVariance">Buffer 'BVariance'.</param>
		/// <param name="result">Modified to contain the outgoing message</param>
		/// <returns><paramref name="result"/></returns>
		/// <remarks><para>
		/// The outgoing message is the factor viewed as a function of 'product' conditioned on the given values.
		/// </para></remarks>
		public static VectorGaussian ProductAverageConditional(Matrix A, [Fresh] Vector BMean, [Fresh] PositiveDefiniteMatrix BVariance, VectorGaussian result)
		{
			// P.mean = A*B.mean
			// P.var = A*B.var*A'
			// if A is invertible, then
			// P.prec = inv(A)'*inv(B.var)*inv(A)
			// P.precTimesMean = inv(A)'*B.precTimesMean
			Vector rmean = A * BMean;
			PositiveDefiniteMatrix rvariance = new PositiveDefiniteMatrix(result.Dimension, result.Dimension);
			Matrix temp = (A*BVariance).Transpose();
			rvariance.SetToProduct(A, temp);
			result.SetMeanAndVariance(rmean, rvariance);
			return result;
		}

		/// <summary>
		/// VMP message to 'mean'
		/// </summary>
		/// <param name="Sample">Constant value for 'sample'.</param>
		/// <param name="Precision">Constant value for 'precision'.</param>
		/// <param name="result">Modified to contain the outgoing message</param>
		/// <returns><paramref name="result"/></returns>
		/// <remarks><para>
		/// The outgoing message is the factor viewed as a function of 'mean' conditioned on the given values.
		/// </para></remarks>
		public static VectorGaussian MeanAverageLogarithm(Vector Sample, PositiveDefiniteMatrix Precision, VectorGaussian result)
		{
			return SampleConditional(Sample, Precision, result);
		}

		/// <summary>
		/// VMP message to 'Sum'
		/// </summary>
		/// <param name="A">Constant value for 'A'.</param>
		/// <param name="B">Incoming message from 'B'. Must be a proper distribution.  If any element is uniform, the result will be uniform.</param>
		/// <param name="MeanOfB">Buffer 'MeanOfB'.</param>
		/// <param name="CovarianceOfB">Buffer 'CovarianceOfB'.</param>
		/// <returns>The outgoing VMP message to the 'Sum' argument</returns>
		/// <remarks><para>
		/// The outgoing message is a distribution matching the moments of 'Sum' as the random arguments are varied.
		/// The formula is <c>proj[sum_(B) p(B) factor(Sum,A,B)]</c>.
		/// </para><para>
		/// Uses John Winn's rule for deterministic factors.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="B"/> is not a proper distribution</exception>
		public static Gaussian SumAverageLogarithm(bool[] A, [SkipIfUniform] VectorGaussian B, Vector MeanOfB, PositiveDefiniteMatrix CovarianceOfB)
		{
			Gaussian result = new Gaussian();
			// p(x|a,b) = N(E[a]'*E[b], E[b]'*var(a)*E[b] + E[a]'*var(b)*E[a] + trace(var(a)*var(b)))
			Vector ma = Vector.FromArray(A.Select(x => x?1.0:0.0).ToArray());
			// Uses John Winn's rule for deterministic factors.
			// Strict variational inference would set the variance to 0.
			result.SetMeanAndVariance(ma.Inner(MeanOfB), CovarianceOfB.QuadraticForm(ma));
			return result;
		}

		/// <summary>
		/// EP message to 'A'
		/// </summary>
		/// <param name="matrixMultiply">Incoming message from 'matrixMultiply'. Must be a proper distribution.  If any element is uniform, the result will be uniform.</param>
		/// <param name="A">Incoming message from 'A'. Must be a proper distribution.  If any element is uniform, the result will be uniform.</param>
		/// <param name="B">Constant value for 'B'.</param>
		/// <param name="result">Modified to contain the outgoing message</param>
		/// <returns><paramref name="result"/></returns>
		/// <remarks><para>
		/// The outgoing message is a distribution matching the moments of 'A' as the random arguments are varied.
		/// The formula is <c>proj[p(A) sum_(matrixMultiply) p(matrixMultiply) factor(matrixMultiply,A,B)]/p(A)</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="matrixMultiply"/> is not a proper distribution</exception>
		/// <exception cref="ImproperMessageException"><paramref name="A"/> is not a proper distribution</exception>
		public static GaussianArray2D AAverageConditional([SkipIfUniform] GaussianArray2D matrixMultiply, [SkipIfUniform] GaussianArray2D A, double[,] B, GaussianArray2D result)
		{
			int rows = matrixMultiply.GetLength(0);
			int cols = matrixMultiply.GetLength(1);
			int inner = B.GetLength(0);
			if (result == null) result = new GaussianArray2D(rows, inner);
			// sum_{i,j} (m[i,j] - a[i,:]*b[:,j])^2/v[i,j] = 
			// sum_{i,j} (m[i,j]^2 - 2m[i,j]a[i,:]*b[:,j] + a[i,:]*(b[:,j] b[:,j]')*a[i,:]')/v[i,j]
			// meanTimesPrec(a[i,:]) = sum_j (m[i,j]/v[i,j]) b[:,j]
			// prec(a[i,:]) = sum_j b[:,j]*b[:,j]'/v[i,j]
			Vector bj = Vector.Zero(inner);
			Vector mean = Vector.Zero(inner);
			VectorGaussian ai = new VectorGaussian(inner);
			PositiveDefiniteMatrix variance = new PositiveDefiniteMatrix(inner, inner);
			for (int i = 0; i < rows; i++) {
				ai.Precision.SetAllElementsTo(0.0);
				ai.MeanTimesPrecision.SetAllElementsTo(0.0);
				// we are projecting from family of full covariance Gaussians to diagonal
				// covariance, so we should include the context
				for (int c = 0; c < inner; c++) {
					ai.Precision[c, c] = A[i, c].Precision;
					ai.MeanTimesPrecision[c] = A[i, c].MeanTimesPrecision;
				}
				for (int j = 0; j < cols; j++) {
					Gaussian xij = matrixMultiply[i, j];
					for (int k = 0; k < inner; k++) {
						bj[k] = B[k, j];
					}
					if (xij.IsPointMass) throw new NotImplementedException(LowRankNotSupportedMessage);
					ai.Precision.SetToSumWithOuter(ai.Precision, xij.Precision, bj, bj);
					ai.MeanTimesPrecision.SetToSum(1.0, ai.MeanTimesPrecision, xij.MeanTimesPrecision, bj);
				}
				ai.GetMeanAndVariance(mean, variance);
				for (int k = 0; k < inner; k++) {
					Gaussian rik = result[i, k];
					rik.SetMeanAndVariance(mean[k], variance[k, k]);
					result[i, k] = rik / A[i, k];
				}
			}
			return result;
		}

		/// <summary>
		/// Update the buffer 'CovarianceOfB'
		/// </summary>
		/// <param name="B">Incoming message from 'B'. Must be a proper distribution.  If any element is uniform, the result will be uniform.</param>
		/// <returns>New value of buffer 'CovarianceOfB'</returns>
		/// <exception cref="ImproperMessageException"><paramref name="B"/> is not a proper distribution</exception>
		public static PositiveDefiniteMatrix Ebbt(PositiveDefiniteMatrix CovarianceOfB, Vector MeanOfB, PositiveDefiniteMatrix result)
		{
			result.SetToSumWithOuter(CovarianceOfB, 1, MeanOfB, MeanOfB);
			return result;
		}

		/// <summary>
		/// Update the buffer 'MeanOfB'
		/// </summary>
		/// <param name="B">Incoming message from 'B'. Must be a proper distribution.  If any element is uniform, the result will be uniform.</param>
		/// <param name="CovarianceOfB">Buffer 'CovarianceOfB'.</param>
		/// <returns>New value of buffer 'MeanOfB'</returns>
		/// <remarks><para>
		/// 
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="B"/> is not a proper distribution</exception>
		public static Vector MeanOfB([Proper] VectorGaussian B, PositiveDefiniteMatrix CovarianceOfB)
		{
			return CovarianceOfB * B.MeanTimesPrecision;
		}

		/// <summary>
		/// VMP message to 'A'
		/// </summary>
		/// <param name="X">Incoming message from 'X'. Must be a proper distribution.  If uniform, the result will be uniform.</param>
		/// <param name="A">Incoming message from 'A'.</param>
		/// <param name="B">Incoming message from 'B'. Must be a proper distribution.  If any element is uniform, the result will be uniform.</param>
		/// <param name="MeanOfB">Buffer 'MeanOfB'.</param>
		/// <param name="CovarianceOfB">Buffer 'CovarianceOfB'.</param>
		/// <param name="result">Modified to contain the outgoing message</param>
		/// <returns><paramref name="result"/></returns>
		/// <remarks><para>
		/// The outgoing message is the exponential of the average log-factor value, where the average is over all arguments except 'A'.
		/// Because the factor is deterministic, 'X' is integrated out before taking the logarithm.
		/// The formula is <c>exp(sum_(B) p(B) log(sum_X p(X) factor(X,A,B)))</c>.
		/// </para></remarks>
		/// <exception cref="ImproperMessageException"><paramref name="X"/> is not a proper distribution</exception>
		/// <exception cref="ImproperMessageException"><paramref name="B"/> is not a proper distribution</exception>
		public static GaussianArray AAverageLogarithm([SkipIfUniform] Gaussian X, GaussianArray A, [SkipIfUniform] VectorGaussian B, Vector MeanOfB, PositiveDefiniteMatrix CovarianceOfB, /*PositiveDefiniteMatrix Ebbt,*/ GaussianArray result)
		{
			int inner = MeanOfB.Count;
			if (result == null) result = new GaussianArray(inner);
			// E[log N(x[i,j]; a[i,:]*b[:,j], 0)] = -0.5 E[(x[i,j]- sum_k a[i,k]*b[k,j])^2]/0 
			// = -0.5 (E[x[i,j]^2] - 2 E[x[i,j]] a[i,k] E[b[k,j]] + a[i,k] a[i,k2] E(b[k,j] b[k2,j]))/0
			// a[i,k] * (-2 E[x[i,j]] E[b[k,j]] + sum_{k2 not k} E[a[i,k2]] E(b[k,j] b[k2,j]))
			// a[i,k]^2 * E(b[k,j]^2)
			// message to a[i,k] = N(a; inv(prec[i,k])*(sum_j E[b[k,j]]*res[i,j,k]/var(x[i,j])), inv(prec[i,k]))
			// where res[i,j,k] = E[x[i,j]] - sum_{k2 not k} E[a[i,k2]] E[b[k2,j]]
			// prec[i,k] = sum_j E(b[k,j]^2)/var(x[i,j])
			// result.Precision = prec[i,k]
			// result.MeanTimesPrecision = sum_j E[b[k,j]]*res[i,j,k]/var(x[i,j]) 
			//                           = sum_j E[b[k,j]]*(X.MeanTimesPrecision - X.precision*(sum_{k2 not k}))

			var Ebbt = new PositiveDefiniteMatrix(inner, inner);
			// should we be caching this too? 
			Ebbt.SetToSumWithOuter(CovarianceOfB, 1, MeanOfB, MeanOfB);
			//var ma = A.Select(z => z.GetMean()).ToArray(); 
			var ma = new double[inner];
			for (int k = 0; k < inner; k++)
				ma[k] = A[k].GetMean();
			for (int k = 0; k < inner; k++) {
				double prec = 0.0;
				double pm = 0.0;
				prec += Ebbt[k, k] * X.Precision;
				double sum = 0.0;
				for (int r = 0; r < inner; r++)
					if (r != k)
						sum += Ebbt[r, k] * ma[r];
				pm += MeanOfB[k] * X.MeanTimesPrecision - X.Precision * sum;

				Gaussian rk = result[k];
				rk.Precision = prec;
				if (prec < 0)
					throw new ApplicationException("improper message");
				rk.MeanTimesPrecision = pm;
				result[k] = rk;
			}

			return result;
		}