I liked @ttnphns's suggestion (that he made in the comments above) so much, that could not resist from trying it out.
As @ttnphns said, LDA is equivalent to canonical correlation analysis (CCA) between your multivariate data and the set of dummy variables coding class labels. In case of only two groups, there is only one non-redundant dummy variable that takes e.g. value of $0$ for one class and value of $1$ for the second class. Running CCA between this dummy variable and the multivariate data will yield the canonical axis (maximizing correlation with the dummy variable) identical to the LDA axis (maximizing between-class/within-class variance ratio).
If we now replace this dummy variable by your class probability variable that takes values from $[0, 1]$, we can still run CCA analysis and use the first canonical axis. It will not be equal to any LDA axis anymore, because now that there are no clear groups LDA does not make sense.
Here is a figure illustrating it:
Blue dots are one class, red dots are another class, size of the dots represents the probabilistic label (smallest dots for $p\approx 0.5$, largest dots for $p\approx 0$ and $p\approx 1$). Solid line is the LDA axis. Dashed line is the CCA axis, computed using probabilistic labels. It goes diagonally because if you look very carefully you will notice that the probabilistic labels were assigned such that they increase along the diagonal.
Note that in this very simple case of only two classes and consequently only a one-dimensional label variable, CCA and LDA are equivalent to good old multiple regression (regressing label variable $l\in \mathbb R$ on $\mathbf x \in \mathbb R^2$ results in a two-dimensional $\beta$ coefficient that defines an axis in $\mathbb R^2$). I demonstrate this equivalence numerically in my code below.
MATLAB code
%// generate data
n = 100;
x = randn([n*2 2]);
x(1:n, :) = bsxfun(@times, x(1:n, :), [1 4]);
x(1:n, :) = bsxfun(@plus, x(1:n, :), [1 1]);
x(n+1:end, :) = bsxfun(@times, x(n+1:end, :), [1 4]);
x(n+1:end, :) = bsxfun(@plus, x(n+1:end, :), [10 10]);
x = bsxfun(@plus, x, -mean(x));
labels = [ones(n,1); 2*ones(n,1)]; %// group labels (ones and twos)
prob_labels = (labels-mean(labels) * x * [1 1]'); %// probabilistic labels (from 0 to 1)
prob_labels = prob_labels / max(abs(prob_labels)) / 2 + 0.5;
%// LDA via within- and between-class covariance matrices (my own function)
axisLDA = mylda(x, labels);
%// CCA with probabilistic labels
[axisCCAprob, ~] = canoncorr(x, prob_labels);
axisCCAprob = axisCCAprob / norm(axisCCAprob);
%// plot
axesScale = 15;
dotScale = 200;
figure
hold on
axis([-20 20 -20 20])
axis square
scatter(x(1:n,1), x(1:n,2), abs(prob_labels(1:n)-0.5)*dotScale, 'b')
scatter(x(n+1:end,1), x(n+1:end,2), abs(prob_labels(n+1:end,:)-0.5)*dotScale, 'r')
plot(axisLDA(1)*axesScale*[-1 1], axisLDA(2)*axesScale*[-1 1], 'k')
plot(axisCCAprob(1)*axesScale*[-1 1], axisCCAprob(2)*axesScale*[-1 1], 'k--')
We can now check that different approaches mentioned above give identical results:
%// check that different approaches give identical results
display(['LDA, take 1: ' num2str(axisLDA')])
[axisCCA, ~] = canoncorr(x, labels); %// LDA via canonical correlation with dummies
axisCCA = axisCCA / norm(axisCCA);
display(['LDA, take 2: ' num2str(axisCCA')])
axisRegression = regress(labels, x); %// LDA via regression
axisRegression = axisRegression / norm(axisRegression);
display(['LDA, take 3: ' num2str(axisRegression')])
display(['CCA with probabilistic labels, take 1: ' num2str(axisCCAprob')])
axisRegressionProb = regress(prob_labels, x); %// CCA via regression
axisRegressionProb = axisRegressionProb / norm(axisRegressionProb);
display(['CCA with probabilistic labels, take 2: ' num2str(axisRegressionProb')])
Running this code produces the following output:
LDA, take 1: 0.99735 0.072685
LDA, take 2: 0.99735 0.072685
LDA, take 3: 0.99735 0.072685
CCA with probabilistic labels, take 1: 0.68477 0.72876
CCA with probabilistic labels, take 2: -0.68477 -0.72876